Classification of dark pion multiplets as dark matter candidates and collider phenomenology

New confining sectors can contain a set of pseudo-Goldstone mesons that exhibit a complicated structure in terms of stability and relative masses. Stable ones can act as dark matter candidates, while their interactions with the unstable ones determine their relic abundances. The overall structure, by specifying which channels are kinematically forbidden or not, affects the cosmology, constraints and collider phenomenology. In this paper, we present a classification of these pseudo-Goldstone meson structures. We find that the structures can be classified into three categories, corresponding to strong, suppressed and essentially non-existent constraints from indirect detection. Limits on decay lengths of the unstable mesons and dark jet properties are presented for several benchmark models.


Introduction
There is little doubt about the existence of Dark Matter (DM) at this point. A plethora of observations such as galactic rotational curves, large and small scale structures and anisotropies in the cosmic microwave background all point toward its existence. Beside some basic facts such that it must be long-lived, non-baryonic and have an abundance roughly five times that of baryonic matter [1], the exact nature of dark matter is however still unknown. One fairly understudied possibility is that dark matter consists of the stable pseudo-Goldstone mesons of a new confining sector. In such a scenario, a set of fermions charged under a new confining group is introduced. Depending on their quantum numbers, masses and other interactions, some approximate chiral symmetry can exist between them. If the fermions are light enough compared to the confinement scale of the new group, this approximate chiral symmetry will however be broken spontaneously to a smaller subgroup by a condensate of the dark fermions. To each of the spontaneously broken symmetry corresponds a pseudo-Goldstone boson or so-called dark pion. If enough symmetries are left unbroken by the condensate and explicit symmetry breaking, some of the dark pions can be stable and act as dark matter candidates. The other dark pions eventually decay away, but their interactions with the stable ones can potentially explain the dark matter abundance we observe today.
Such a scenario is interesting to study for several reasons. First, dark pions are present in many existing models of beyond the Standard Model physics and are theoretically wellmotivated. They are ubiquitous in models of Neutral Naturalness, of which Twin Higgs [2,3] is arguably the archetype. They also have been proposed as a potential explanation for dark matter [4,5] and certain excesses [6], as well as in more esoteric settings [7]. Second, multiplets of dark pions allow for explanations of the dark matter abundance that fall outside the standard paradigm of thermal freeze out of Weakly Interacting Massive Particles (WIMP) [8][9][10]. Strongly Interacting Massive Particles (SIMP) [11], in which 3 → 2 processes determine the DM relic abundance, is a well-known example of this, but dark pions also naturally accommodate more exotic mechanisms such as Codecaying Dark Matter [12][13][14] and ELDER Dark Matter [15]. Such scenarios are typically far less constrained by direct detection than WIMPs. Third, there has recently been an increasing amount of work dedicated to discovering jets of dark hadrons, or so called dark jets. From a theoretical point of view, this includes emerging jets [16], semivisible jets [17] and softbombs [18]. From an experimental point of view, collider searches have been performed in ref. [19]. Of course, one of the practical problems of looking for dark jets is that they can be realized at colliders in many different ways, depending on their multiplicity, thrust and the stability of their constituent dark hadrons. Being able to narrow down what dark jets should look like would certainly simplify experimental search strategies.
With this context in mind, the goal of this article is to classify multiplets of dark pions in terms of their symmetry breaking structure, analyze their resulting viability as dark matter candidates and map this to the corresponding dark jet properties at colliders. More precisely, the symmetry breaking takes place in two steps: Spontaneous Chiral Symmetry -2 -

JHEP02(2020)196
Breaking (CSB) of a flavor group to a smaller subgroup and Explicit Breaking (EB) of the latter by interactions or masses. The combination of these two breakings results in a set of dark pions with a possibly complicated structure of mass versus stability. This structure crucially affects the evolution of the dark matter abundance in the early Universe. Making sure that the resulting dark matter is compatible with the Universe as we observe it today imposes constraints on the structure of the dark pion sector. This in turn is reflected in how dark jets should look like at colliders. What we do in this article is classify these structures, apply cosmological constraints and make predictions as to the characteristics of dark jets. More specifically, we will concentrate on multiplets that contain both stable and unstable dark mesons at the same time. If all dark pions were stable, this would reduce to SIMP or possibly standard WIMP thermal freeze-out, while, if they were all unstable, they obviously could not serve as dark matter candidates. In a lot of ways, this paper is a generalization of ref. [20]. Previous work on combinations of stable and unstable dark mesons also include refs. [6,14,[21][22][23].
The final result of this paper will be that dark meson multiplets can be separated into three distinct categories. In the first category, all stable pions can annihilate with another to produce unstable ones, even at low temperature. This results in strong bounds from indirect detection and jets that are similar to QCD jets with a higher confinement scale. In the second category, some stable pions cannot annihilate with other ones to produce unstable pions, but some others can. This considerably alleviates the bounds from indirect detection, opening the way to less QCD-like jets. When the mass splitting between the dark pions is small, the dark jets can take many forms. However, a larger splitting between the dark pions forces the jets to be increasingly narrow, with a larger proportion of stable mesons, larger multiplicity and with short decay lengths for their unstable dark pions. In the third category, all annihilations of two stable pions to a final state containing an unstable pion are forbidden. This means that bounds from indirect detection are essentially inexistent. The resulting dark jets are similar to those of the second category, but with a larger upper limit on the decay length of the unstable mesons.
The article is organized as follows. First, we classify the different patterns of symmetry breaking. The mechanism through which they can lead to the correct dark matter relic abundance is then described. The process through which we obtain current cosmological bounds is explained afterward. Limits on the mass splitting and decay lengths are then presented for a set of benchmark models. Collider properties such as the fraction of the dark jets that consists of stable mesons and the thrust and multiplicity of the dark jets are discussed considering these constraints. We finish with some concluding remarks. The first appendix elaborates on the technical details of the procedure through which we obtain the correct dark matter abundance. The second appendix discusses more carefully our benchmark symmetry breaking patterns.

Concepts and classification of symmetry breaking patterns
We begin this article by introducing the concepts and notation that we will use throughout the paper. We will then classify the different combinations of chiral and explicit symmetry breaking.

Notation
First, we define more clearly what we mean by a new confining sector. We begin by assuming a new confining gauge group G. Considering the goals of this paper, we will generally assume its confinement scale is at a value that can reasonably be probed at current colliders, but this is not crucial from a cosmological point of view. We then introduce a set of N massive Dirac dark quarks q i of mass m i . By definition, they are charged under G but neutral under SM gauge groups. We assume there is an at least approximate flavor global symmetry between the dark quarks characterized by the group G. We also assume that the dark quarks communicate with the SM sector, but that the interactions between the two sectors are rather weak.
Next, we consider three patterns of chiral symmetry breaking. They are: The groups on the left G represent the original flavor symmetries and the groups on the right H what they are broken to by the condensate. The first pattern is typical of dark quarks charged under complex representations of G, the second to pseudoreal representations and the third to real representations [24]. To each broken symmetry corresponds a pseudo-Goldstone meson or dark pion. Together, they form a representation of H whose size is included in eq. (2.1). To alleviate the text, we will in general refer to dark pions as simply pions, as there will be no risks of confusing them with the usual SM pions. They are referred to as pseudo-Goldstone bosons as we do assume the presence of mass terms that render the spontaneously broken symmetries approximate and give mass to the pions. The pion decay constant associated to the pions is labelled as f . The advantage of including real and pseudoreal representations for dark quarks is that it opens many possibilities for the structure of the dark sector. Perhaps most interestingly is that dark color singlets can then be constructed with two quarks or two antiquarks, not necessarily a quark and an antiquark. This is specially useful for real dark quarks, where a dark color singlet can be constructed with twice the same quark, resulting in a meson with a mass independent of the masses of the other dark quarks. Then, we assume that the flavor group H that was left unbroken by the condensate is explicitly broken by masses or interactions to some smaller subgroup h. The pions will then decompose themselves into a set of representations of h. This symmetry will generally be sufficient to maintain some of them stable, but not all of them. We will assume that no particles beside the pions exist that are charged under this group. Otherwise, the stability of the pions would not be insured anymore and this would effectively be equivalent to the symmetry being broken. The breaking of H will also in general result in the different multiplets of h having different masses.
To summarize, the pattern of symmetry breaking is:

SU (4 )→SO (4 )→U (1 )
(c) Category III. The first breaking determines the number of pions and the second breaking their masses and stability. The combination of these breakings determines a structure of pions with given masses and stability. This structure in turn controls which processes are kinematically forbidden or not. This determines the cosmology and the collider phenomenology.

Classification of symmetry breaking patterns
With the notation established, the structures of pions resulting from a given combination of chiral and explicit symmetry breaking can be classified into three categories from both a theoretical and phenomenological point of view. In the hope of making the discussion more concrete, we will provide an example for each case.
Category I: unforbidden. A structure of pions is said to belong to category I if, for each stable pion, there exist another stable pion with which it can annihilate to produce at least one unstable pion via a kinematically allowed 2 → 2 pion scattering process. 1 In more casual terms, it simply means that every stable pion can annihilate with another to produce at least an unstable one, even at very low temperature.
An example of such a spectrum is shown in figure 1a. In it, three dark quarks transform under a complex representation of G. This results in an SU(3) × SU(3) → SU(3) pattern of chiral symmetry breaking, leading to eight pions. The resulting SU(3) is then explicitly broken to SU(2) × U(1). This splits the pions into multiplets (3, 0), (2, ±1) and (1, 0). Assuming the dark quark not charged under SU(2) is q 3 and that it is lighter than the two others, the lightest state is then the singlet. Since it does not carry any charge under SU(2) × U(1), this group is insufficient to maintain the singlet stable and it will in general decay, barring any additional symmetry. All other pions can be proven to be stable because of conservation of group charges and energy. By the same principles, it is easy to verify that every stable pion can annihilate with its conjugate to produce a pair of singlets, thereby satisfying the condition of category I.

JHEP02(2020)196
Category II: partially forbidden. A structure of pions is said to belong to category II if, for some but not all of its stable pions, there exist another stable pion with which they can annihilate to produce at least one unstable pion via a kinematically allowed 2 → 2 pion scattering process. Simply put, there are some stable pions that can collide with some other stable pions to produce unstable pions and some that can't.
An example of such a spectrum is shown in figure 1b. It is similar to the example of category I, but with m 3 > m 1 instead. It is simple to use group charges and energy conservation to prove that the triplet and doublets are once again stable. It is also easy to verify that the (3, 0) cannot produce unstable pions via a kinematically allowed annihilation with another stable pion. However, the annihilation of (2, +1) with (2, −1) to (3, 0) and (1, 0) is kinematically allowed and has a non-zero amplitude even at leading order. All in all, this means that only a fraction of the stable pions can produce unstable pions at very low velocity.
Category III: forbidden. A structure of pions is said to belong to category III if, for every stable pion, there does not exist another stable pion with which it can annihilate to produce at least one unstable pion via a kinematically allowed 2 → 2 pion scattering process. In simple terms, it means that it is impossible for two stable pions to produce an unstable one via a collision at low temperature.
An example of such a spectrum is shown in figure 1c. It consists of two quarks transforming under a real representation of G. The pattern of spontaneous symmetry breaking is then SU(4) → SO(4), resulting in nine pions. The group SO(4) is then broken explicitly to its U(1) subgroup in which dark quarks have charge +1/2 and dark antiquarks −1/2, i.e. the equivalent of the baryon number up to a normalization chosen for convenience. Pions made of a quark and an antiquark therefore have a charge of 0, while those made of two quarks or two antiquarks have respectively a charge of +1 or −1. Only the pion made of twice the lightest quark and its conjugate are therefore stable. They also cannot annihilate between each other to produce unstable pions via kinematically allowed processes.
Additional comments. This classification is very pragmatic. Indeed, the annihilation of stable pions to unstable ones provides constraints from indirect detection. For category I, all stable pions can contribute to the indirect detection signal, which results in strong constraints. For category II, only a subset of stable pions can contribute significantly to the indirect detection signal, which significantly reduces the bounds from indirect detection. For category III, none of the stable pions can contribute significantly to the indirect detection signal, which results in essentially non-existent constraints. Of course, these categories also affect the cosmological evolution.
Radiative corrections can potentially determine in which category a set of pions belongs to. For example, radiative corrections make the charged pion heavier than the neutral one in the Standard Model. Ultimately, the only factor that matters is the masses of the pions, not whether they come from the quark masses or radiative corrections. To avoid any potential confusion, we will only use benchmarks for which the category is not determined by radiative corrections given sufficient mass splitting between the dark quarks. As such, JHEP02(2020)196 Table 1. Benchmark scenarios for the two categories of pion structures. The condition represents a relation that is assumed between the dark quark masses. It is sometimes crucial for the structure to be in a given category but not always. See appendix B for more details.
we will also neglect all radiative corrections. Doing so would anyhow require knowledge beyond the simple dark quark content of the new confining sector. In some cases, different multiplets can be degenerate in mass. Unless there is some symmetry involved, radiative corrections will generally lift this degeneracy. Scattering between such particles are typically impeded by a small phase-space.
A series of examples that will also serve as benchmarks is shown in table 1. We refer to appendix B for a complete description of every benchmark model. Note that discrete subgroups could of course also be used, even though we did not include any examples.

Overview of the cosmological evolution
Having established the notation, we now proceed to discuss the cosmological evolution of the pion densities. This section will include a discussion of our assumptions, the different regimes of density evolution and how the pion structures affect the abundances.

Assumptions
We begin this section by stating the assumptions we will make during the computation of the dark matter relic abundance. None of them are crucial to our discussion and are simply meant to simplify the presentation. First, we assume that annihilation of pions directly to SM particles is negligible compared to annihilation to other pions. Otherwise, bounds from direct detection would be difficult to satisfy (see ref. [20]). Dark pions would also act mostly as standard WIMPs, which are already well studied. Second, we assume that the dark sector maintains the same temperature as the SM sector. This simplifies the presentation tremendously as there is a very large number of ways for the two sectors to exchange energy. In practice, requesting that the mesons decay fast enough to SM particles implies a minimum amount of interaction between the two sectors. If the pions interact with -7 -JHEP02(2020)196 light particles, it should be sufficient for the two sectors to maintain kinematic equilibrium (see ref. [20]), but could prove to be insufficient if the pions interact exclusively with heavy SM particles. One way or the other, this should not affect qualitatively any results. Third, we will neglect the presence of other composite particles. As they should be heavier than the pseudo-Goldstone mesons, they should be negligible in most regions of parameter space and not change the conclusions in the small regions where their effects become important. 2

Different evolution regimes
With these assumptions established, the cosmology can fall into two distinct regimes which we present here. This discussion mirrors the previous work of ref. [20]. Examples of the thermal evolution are shown in figure 2 for a benchmark of each category. TheŶ i 's represent the number density of a given multiplet divided by the entropy density and the size of the multiplet. We refer to appendix A for the mathematical details.
Codecaying Dark Matter regime. The first regime is characterized by the unstable pions having small decay widths, where what is meant by small will soon be clear. The evolution of the dark matter goes through (up-to) three steps as it passes from small to large x ≡ m π 0 /T , where π 0 is some arbitrary reference pion. This is illustrated in figures 2a, 2c and 2e.
First, the Wess-Zumino-Witten (WZW) term allows for the presence of 3 → 2 processes [28][29][30]. This term is present as long as π 5 (G/H) = Z. This is satisfied for all patterns of chiral symmetry breaking of eq. (2.1) with at least two dark quarks, except SU(2) × SU(2) → SU(2) which we will not be using anyhow. At small enough x, 3 → 2 processes are sufficient to maintain the pions at their thermal equilibrium densities. This is the basic mechanism behind SIMP dark matter. If the WZW term is not allowed, this step is simply absent.
Second, as x increases and the pion densities decrease, 3 → 2 processes eventually become inefficient and pions become unable to maintain their thermal equilibrium densities. The 2 → 2 processes are still efficient however and the pions therefore maintain chemical equilibrium with each other. The unstable ones decay and are replaced by the annihilation of stable pions to unstable ones. This leads to an overall decrease of the pion abundance.
Finally, the pions decouple from each other. The end result is a net density of stable pions, while the unstable ones simply decay away.
For a fixed mass structure, two sets of parameters affect the relic abundance. First, decreasing the pion decay constant allows the pions to maintain chemical equilibrium for a longer time, allowing more stable pions to be converted to unstable ones which then decay. In addition, it renders the 3 → 2 processes more efficient. Both of these factors lead to a smaller DM relic abundance. Second, larger pion decay widths mean that the unstable pions will decay faster and will again reduce the DM relic abundance.
Coupling-independent regime. As one increases the decay widths of the unstable pions, their densities approach their equilibrium values. Eventually, a point is reached JHEP02(2020)196 (a) CIa: Codecaying. 10 -14 10 -14 10 -14 10 -14 10 -14  where unstable pions are still in thermal equilibrium when the stable pions decouple from them. Past this point is the so-called coupling-independent regime of ref. [20], dubbed this way because further increasing the decay widths of the unstable pions does not affect the relic densities anymore. Effectively, the unstable pions become just another component of the plasma. This regime is illustrated in figures 2b, 2d and 2f. For a fixed mass structure, the only factor that affects the DM relic abundance is how much the pions interact with each other, i.e. the pion decay constant f .

Cosmological evolution for the different pion structures
The evolution of the pion densities follows the order of the previous subsection for all pion mass structures. The final relative abundances of the pions are however strongly dependent on the category.
In category I, the annihilation of two stable pions to at least an unstable one is always kinematically allowed. This means that all stable pions have similar cross sections for production of unstable pions, up to phase-space factors and numerical coefficients of O(1). This signifies that all the stable pions will decouple from the unstable ones at a similar temperature. This results in all the stable pions having similar densities, as can be seen in figures 2a and 2b.
In category II, only some of the stable pions can annihilate to at least an unstable pion via a kinematically allowed process. Those that are not allowed to do so will have cross sections for the annihilation to unstable ones that are Boltzmann suppressed. This results in them decoupling when their densities are much larger than those that can annihilate to unstable ones. The end result is that the stable pions that cannot annihilate to unstable ones dominate the relic density, as can be seen in figures 2c and 2d.
In category III, none of the stable pions can annihilate to at least an unstable pion via a kinematically allowed process. The lightest pions typically dominate the relic density, as can be seen in figures 2e and 2f.

Constraints
In this section, we discuss the procedures through which constraints are applied and we refer to section 5 for limits on different benchmark models. In this work, we use the experimental value for the relic abundance Ω obs h 2 = 0.1200 ± 0.0012 at 68% confidence [1].

Indirect detection
In models where the dark matter populates rich hidden sectors, annihilations may proceed via sequential decays, giving rise to multibody final states. In fact, depending on the spectrum of the dark sector, stable dark particles colliding with each other may produce unstable dark states. These subsequently decay back into Standard Model particles, which may in turn decay until only stable particles like positrons, antiprotons and photons are left. These particles can potentially be detected and/or affect their environment. Their indirect detection signals depend essentially on the annihilation rate of particles within the hidden sector.

JHEP02(2020)196
The injection of photons and other high energy secondary particles is constrained by the measurement of the CMB by Planck [31], the bounds set by the Fermi-LAT collaboration from DM searches in the Dwarf Spheroidal Galaxies of the Milky Way [32] and by measurements of the positron flux by AMS-02 [33,34]. As a consequence, constraints from the above experiments may provide an upper bound on the self-interaction strength of the pions. These cascade decays have been studied in ref. [35]. We compute our limits using their bounds on bb final states, which should be rather accurate for all hadronic decays and correspond to 95% confidence. Our constraints will be computed using an effective cross section defined by where Y tot is the total number density of pions per entropy density, Ω tot the total relic density, Ω obs its observed value, P 2→2 the set of all distinct 2 → 2 scattering between representations,Ŷ A the number density of pions from a representation A per entropy density divided by the size of the representation and #π P unstable out the number of unstable pions in the out state of a process P . See appendix A for more details. It is obtained by simply combining the rates of unstable pion production from all channels. This expression assumes the proportionality of local and global relative abundances, or so-called proportionality ansatz [36][37][38][39][40][41], which should work to great approximation in the mass range we study. It is also possible to set limits studying the indirect detection signals from the annihilation of two dark pions directly to SM particles, although they are typically weaker than bounds from direct detection if we assume rather weak interactions between the SM and the dark sector [20].
Finally, it is expected that the indirect detection bounds will change in the near future. The lack of a signal from Planck, Fermi-LAT and AMS-02 will lead to tighter limits, but not change the results qualitatively [35]. An improvement of roughly one order of magnitude for the Fermi-LAT bounds will be possible depending on new dwarf spheroidal galaxy discovery [42]. Finally, future more powerful instruments, such as GAMMA-400 [43,44], HERD (High Energy cosmic Radiation Detection) [45,46] or CTA (Cherenkov Telescope Array) [47][48][49] may set even stronger limits at different mass scales.

Other constraints
In our convention, an approximate perturbative unitarity limit of m π /f 4π can be set. As a consequence, as we will see in section 5, there will be an upper bound on the decay lengths of unstable pions. This bound will be the stronger limit in benchmarks where the effective thermal averaged cross section is suppressed or vanishing (category II and III). However, this bound should be taken more as a rough estimate, as chiral perturbation theory breaks down around this limit and additional resonances must be taken into account.
Furthermore, the decay of the unstable pions may disturb Big Bang Nucleosynthesis (BBN). In fact, if the unstable pions decay hadronically, they will inject relatively longlived hadrons that can modify the ratio of protons to neutrons. As a consequence we must require an upper bound on the decay time of about 0.1 second for a 95% confidence -11 -
Finally, the assumption that the annihilation of pions directly to SM particles is negligible compared to the annihilations between pions makes direct detection constraints irrelevant.

Limits on mass splitting and decay lengths
In this section, we present limits on decay lengths as a function of the mass splitting for all three possible structure categories. Of course, it is impossible to consider every pattern of symmetry breaking for every category. We will therefore present a series of benchmarks that illustrate general features. These features should apply roughly to all models within a given category.

Limits on decay lengths for category I
Category I is characterized by all stable pions being able to contribute to the indirect detection signal. This simplifies the analysis, but also means that bounds from indirect detection are very strong.
One unstable pion. We begin by considering patterns with only one unstable pion. The pattern CIa is a good example of this. This benchmark is characterized by five parameters: where m π s 3 is the mass of the triplet and Γ π d 1 the decay width of the unstable pion. We will fix f by requiring that the abundance of the stable dark pions corresponds to the currently observed dark matter abundance. We set N c to 3, as we will do for all SU(N ) cases in this paper. The number of colors N c will be set to 4 for all Sp(2N ) and SO(2N ) cases. 3 This parameter is only relevant for the 3 → 2 processes, which generally do not change qualitatively the final abundances. We will also swap the decay width for the decay length for practicality. This conveniently leaves us with only three parameters.
We show in figure 3 the allowed parameter space for this benchmark as a function of m 3 /m 1 and the decay length of π d 1 for different fixed m π s 3 . The figure also contains contours of constant m π s 3 /f . At small decay length, these are vertical as the system is in the coupling-independent regime. This ratio starts to increase abruptly once the system enters the coupling independent regime, as f must be made smaller to compensate for a too small decay width. Eventually, the pions are required to interact too much with each other and the system becomes incompatible with the limits from indirect detection. As can be seen, the region of m 3 /m 1 close to 1 is always unexcluded. This is simply because in this limit the s-component of the effective cross section goes to zero. As the mass of the triplet increases, the allowed region of parameter space increases. Do note that a large part of the coupling-independent region of parameter space is very close to the current experimental limit from indirect detection. To illustrate this, we include blue dashed curves that indicate what the limits from indirect detection would be if the limit on the effective cross section were to change by ±10%. Obviously, improving the bounds from indirect detection would massively reduce the amount of parameter space available. One can see that the ratio m Two unstable pions. We then consider patterns with two unstable pions. The pattern CIc is a good example of this. This benchmark is characterized by seven parameters (see appendix B for notation): 2) The parameters f and N c are set as before, which leaves us with five parameters. For illustration purposes, we set m 2 /m 1 to 1 and assume that Γ π m 1 and Γ π m 2 are proportional, leaving three parameters. The results are shown in figure 5. At m 3 /m 1 close to 1, the upper limit on the decay length of the lightest unstable pion is raised in comparison to benchmark CIa. This is simply because the presence of another unstable pion means that the lightest one doesn't need to decay as fast. As m 3 /m 1 approaches 0, the difference between these two cases vanishes as the heaviest unstable pion is simply not present in sufficient quantity to affect the relic abundances.

Limits on decay lengths for category II
Category II is characterized by only a subset of the stable pions being able to contribute to indirect detection. This generally results in suppressed bounds from indirect detection.
One unstable pion. We begin this section by considering examples with only one unstable pion. Benchmark CIIb is a good example of this. As this is simply benchmark CIa with the mass hierarchy inverted, the same discussion about the parameter space applies. The constraints are shown in figure 6. As can be seen, a larger decay width or a larger mass splitting must be compensated by a smaller pion decay constant, which eventually becomes incompatible with unitarity. The excluded corner in the upper left comes from indirect detection and the slowly dropping curve on the right from unitarity. The indirect detection signal is suppressed in most of the parameter space because a sizable mass splitting reduces the abundance of the doublets which in turn reduces the effective indirect detection cross section. There is an upper limit on the decay length of the unstable pion. For the range of pion masses we consider, it is of O(10) to O(100) m and decreases as the pion masses increase. However, such an upper limit is only possible close to a very specific value of m 3 /m 1 , which seemingly does not have any special theoretical justification. In contrast to category I, a much larger range for the pion decay constant is possible. Between      Two unstable pions. A good example of a pion structure of category II with two unstable pions is the benchmark CIIe. This structure is characterized by 7 parameters: 3) The parameters f and N c are set as before. For the sake of the presentation, we assume m 3 m 1 and m 4 m 1 to be equal and Γ π m 1 and Γ π m 2 to be proportional. The results are shown in figure 8. As can be seen, the presence of an extra unstable pion makes a small difference at small m 3 m 1 , but this effect quickly disappears because of Boltzmann suppression.

Limits on decay lengths for category III
Finally, structures of category III are those for which all annihilations of two stable pions to at least an unstable one are kinematically forbidden. These are very similar to category II, except for the fact that the bounds from indirect detection essentially do not exist anymore. As such, we will present only one example, which we take to be CIIIa. It is characterized by 11 parameters:       and seven decay widths. We set f and N c as usual. For the sake of presentation, we assume the decay widths to be proportional to each others, with the decay width of π d 1 to be the largest. The results are shown in figure 9. As expected, they are similar to those of category II without the bounds from indirect detection. This allows for considerably longer decay lengths.

Dark jet properties
In this section, we investigate the properties of dark jets. In particular, we present results for the jet multiplicity, the fraction of unstable particles in a dark jet (r inv ) and the thrust. We qualitatively outline how these variables vary due to the different features of the models within a given category. In the following, we fix the value of the confinement scale as

Multiplicity
The parton shower is followed by a non-perturbative process, the fragmentation process, and there is a relation between the number of the radiated partons and the number of the produced hadrons. The average number of hadrons is just the first Mellin moment of the fragmentation function [59]. In the modified leading-log approximation, the hadron -19 -
Excluded by Perturbative unitarity Figure 9. Allowed region of parameter space for the benchmark CIIIa for different ratios of the decay widths. The mass of the lightest pion is kept at a constant value of 75 GeV. The contour lines represent constant values of m π a 1 /f . The red region is excluded by perturbative unitarity constraints. The decay length is that of the π d 1 .
multiplicity is calculable up to a normalization constant [59,60] where A is the normalization constant, C F and C A are the Casimir invariants for the fundamental and adjoint representations, while T F is the normalization of the generators Tr(t a t b ) = T F δ ab . We normalize the generators of all the groups in the standard way T F = 1/2. Therefore, for the different gauge groups, we have Furthermore, b 1 is the one loop coefficient of the beta function, and α s (ŝ) = (b 1 logŝ/Λ 2 d ) −1 is the running strong coupling constant for a centre of mass energy squaredŝ. The one loop coefficients of the beta function are: (6.4) This prediction for the average multiplicity is valid for sufficiently large values of the coefficient b 1 and has been verified experimentally for QCD [61]. Notice also that the group O(N c ) stops confining for a smaller number of fermions with respect to SU(N c ) or -20 -

JHEP02(2020)196
Sp(N c ). Furthermore, its one loop coefficient of the β function is smaller than the one of the other groups and therefore it radiates more.
The average multiplicity depends on the number of colors N c , the number of flavors N and the ratio between the scale of the centre of mass energy and the confinement scale √ŝ /Λ d . For larger N , the running of the dark sector coupling to smaller values is slower and more partons are radiated at higher scales. This results in a larger number of dark pions. In figure 10, we show the behavior of the average multiplicity as a function of the ratio √ŝ /Λ d for SU(3) (10a), Sp(4) (10b) and O(4) (10c). Notice that we keep the same normalization factor for each curve, A = 0.1. The normalization factor A cannot be theoretically determined unless we know how gluons fragment into hadrons. Since it can only be obtained by fit to experimental data, 4 these curves should be taken as qualitative behaviors rather than explicit values. Increasing the number of dark colors N c has similar effects on the average multiplicity to lowering N , since they both enter the one loop coefficient of the beta function b 1 . Figure 10 shows that the average multiplicity is also strongly dependent on the value of the dark confinement scale Λ d , that enters logarithmically in the strong coupling constant α s (ŝ). In particular, the larger is Λ d , the smaller is the average multiplicity of dark mesons in a jet produced in a collision with centre of mass energy √ŝ . If we assume that Λ d ∼ f / √ N c , we expect to have larger multiplicity for smaller f . This means that for the models of category I with one unstable pion, we have larger multiplicity for smaller dark quark splitting (compare figures 3 and 4). Notice that this conclusion changes if we have more unstable pions, depending on the parameters of the model, see figure 5. On the other hand, for category II and III the multiplicity increases for larger mass splittings (see figures 6, 7, 8 and 9).

Average collinear missing transverse energy fraction, r inv
Generically, both stable and unstable dark pions can be produced in a parton shower. Assuming the unstable pions decay hadronically, this will result in the production of semivisible jets [17,62,63]. These roughly look like generic jets, but are accompanied by collinear stable particles that escape the detector and contribute to the missing transverse energy (MET). A useful quantity is therefore r inv , defined as the average fraction of the energy of a given semivisible jet that is transmitted to MET. This variable can take on any value between zero (no MET present) and one (no visible jets present). The purpose of this section is to provide a mapping between the allowed parameter space obtained by cosmological observations and the variable r inv .
In the following, we will neglect the contribution of dark baryons. This is justified by the fact that in an Sp(N c ) theory there are no stable baryons, since a baryon would decay into N c /2 mesons [30]. In theories where stable baryons do exist, their production is suppressed at the 10% level in QCD and should be further suppressed in the large N c limit [64]. On the other hand, vector meson production might not be overly suppressed. Their stability will however not in general be insured by any symmetry. Under our assumptions, they should mostly decay quickly to dark pions, in a similar fashion to how SM 4 For example, for QCD with ΛQCD = 226 MeV we would have A 0.12 [61].  vector mesons decay mainly to SM pseudoscalar mesons, which will generally be kinematically allowed. 5 As a consequence, the quantity r inv is roughly controlled by the fraction of stable pions modified by suppression factors for the production of heavier mesons.
In the Lund string model for fragmentation, the production of heavier mesons during hadronization is exponentially suppressed [68]. This is described by the suppression factor for estimating the ratio of q i to q j production where m i > m j are the masses of the dark quarks and Λ d is the dark confinement scale. The larger the mass splitting between q i and q j compared to the confinement scale Λ d , the more suppressed will be the production of dark mesons containing the heavier quark. As a consequence, this may reduce the number of stable states in the dark pion shower. Notice that the values of the quark masses are related to m π and B 0 . Consider the spectrum of the case CIa (figure 11a). In this scenario, we have m 1 = m 2 = m 2 π s 3 /(2B 0 ). Therefore, the suppression factors are where r = m 3 /m 1 . For a small mass difference between the dark quarks m 1 = m 2 and m 3 , the factor T 13 in eq. (6.6) will be close to one and therefore we expect equal production of all three dark quarks. The value of r inv will be close to the ratio between the number of stable and unstable states, given roughly by r inv ∼ 7/8 (see tables in appendix B.2 for details on the stability structure). As the mass splitting increases, the production of the mesons containing only q 1 or q 2 is suppressed with respect to the mesons containing q 3 . As a consequence, we expect that T 13 decreases for larger mass splitting. How fast T 13 decreases depends on the value of the quark masses, or equivalently the ratio between m π and B 0 . Figure 11a shows that for large splitting and B 0 = m 1 we obtain T 13 < 0.1. On the other hand, a different value of B 0 would produce slightly different results. Choosing the more sensible B 0 = f , the suppression factor T 13 decreases slowly reaching a value 0.8 for large mass splitting. This happens because the increasing value in the mass splitting (r 1) is almost compensated by the increasing value of f . This is shown in figure 11b. In either cases, a larger mass splitting implies more visible jets.
Interestingly, case CIIb (with the same breaking pattern, but with a different mass hierarchy m 1 = m 2 < m 3 ) leads to a different behavior. For small mass splitting, T 31 ∼ 1 and r inv ∼ 7/8 as for case Ia. However, as the mass splitting increases, T 31 tends to 0 very quickly, leading to very invisible dark jets. The behavior is similar for both B 0 = m 1 and B 0 = f . All models of category I are expected to have a similar behavior as that of benchmark CIa. Conversely, all models of category II and III have similar behaviors to those of JHEP02(2020)196 benchmark CIIb. The exact range over which r inv can vary however depends on the details of the model.

Thrust
Shape variables [59] are a common approach in order to study the jet-like properties of hadronic final states. These quantities can characterize whether the distribution of the hadrons in a jet is pencil-like, planar or spherical. This is a strong way to obtain information on a jet. A widely used quantity is the thrust [59,69], defined as where p i are the final-state hadron momenta and n is a unit vector. This variable is bounded between 1/2 ≤ T ≤ 1, where a pencil-like event has T = 1, while a spherical one has T = 1/2. The thrust is infrared and collinear safe: an additional parton that is soft or radiated collinear to the thrust axis will not change the thrust. Its distribution is well described by [70] 1 where √ŝ is a high-energy scale associated with the total system being probed by the thrust and C F is the Casimir invariant. We listed the Casimir invariant for the different gauge groups in eq. (6.3).

JHEP02(2020)196
The thrust depends on the number of colors N c , the number of flavors N , the confining scale Λ d and the centre of mass energy √ŝ . In figure 12, we show 1 − T as a function of the ratio √ŝ /Λ d , for different values of N and for the different groups SU(3) (a), Sp(4) (b) and O(4) (c). The thrust for the groups SU(3) and Sp(4) gives similar results. This is due to the fact that the running of their gauge coupling α s is similar. In particular, SU (3) with N = 3 and Sp(4) with N = 3 have the same coefficient b 1 . On the other hand, models with a gauge group O(N c ) have a larger strong coupling. This results in generally more spherical jets for O(N c ) models with respect to SU(N c ) or Sp(N c ) models. Furthermore, the figure shows that a larger confining scale Λ d would lead to smaller thrust and therefore more spherical objects.

Conclusion
In this work, we investigated models where the dark matter consists of the stable pseudo-Goldstone mesons of a new confining sector. The combination of spontaneous and explicit symmetry breaking results in a complicated structure of pions in terms of stability and masses. This structure determines which processes are kinematically forbidden or not. In turn, this affects the cosmological evolution, the current constraints and the range of possible collider signatures. We found that all pseudo-Goldstone meson structures with both stable and unstable pions can be classified into three categories.
The first category includes models where all stable dark pions can annihilate with other stable pions to produce unstable ones. As a consequence, all stable pions have similar cross section for the production of unstable ones and therefore decouple from the unstable pions at a similar temperature. As a result, they all have similar densities and there are strong bounds from indirect detection dark matter searches. The jet properties of dark jets cannot vary over a large range and are not that different from QCD jets with a higher confinement scale. At least one unstable pion is required to have a decay length of less than O(1) m for the masses considered in our benchmarks. These results lead to many different possible signatures at colliders, such as emerging or semivisible jets.
In the second category, not all the stable pions can annihilate with other stable pions to produce unstable ones. Those that are not allowed to do so will have Boltzmann suppressed cross sections for the annihilation to unstable pions, resulting in an earlier decoupling. The dark matter abundance will then be dominated by those pions. Since they do not contribute to the indirect detection signal, the bounds from indirect detection searches are considerably suppressed. The shape of the dark jets, on the other hand, can take many forms, depending on the dark quark mass splitting. Larger splittings will in general require smaller values of f in order to obtain the correct dark matter abundance (and therefore a smaller confining scale Λ d ). As a result, dark jets will be increasingly narrow, with more stable mesons, larger multiplicity and shorter decay lengths. For smaller splittings, considerably longer decay lengths are allowed for the unstable pions, with some benchmarks allowing for decay lengths of up to O(100) m.
Finally, the third category includes models where all the annihilations of two stable pions to unstable ones are forbidden, resulting in the relic density being dominated by the

A Relic dark matter abundance calculation
We present in this appendix all the information relevant for computing the dark matter relic abundance in the examples considered in this paper.

A.1 Lagrangian
The leading order chiral Lagrangian is given by: with Π = π a t a . The matrices t a correspond to the 2N × 2N broken generators (see section B.1 for their expressions and normalization). The B 0 parameter is simply a constant that can be traded in practice for the mass of one of the pions. In all cases we consider in this paper, M can be written as: where M is a diagonal matrix whose entries satisfy M ii = m q i . In addition to L 0 , the Wess-Zumino-Witten term can possibly be present:

JHEP02(2020)196
A.2 2 → 2 scattering Consider a 2 → 2 pion scattering process p given by p : π a π b → π c π d . (A.5) The amplitude of p is with the coefficients where the sums are over all permutations of the indices inside the brackets, even repeated ones. For example, eq. (A.6) can be rewritten as: (A.12) One can then define an effective amplitude square as The first possible factor of 1/2 in eq. (A.14) is to avoid double counting collisions between identical particles. The second factor is to take into account identical particles in the final state such that all integrals over the polar angle can be performed from 0 to π.
A great way to simplify the numerical computations is to use the fact that pions inside an unbroken multiplet maintain the same density at all times. Consider four nonnecessarily distinct multiplets A, B, C and D. Define P : π A π B → π C π D as the set of distinct p : π a π b → π c π d such that π a ∈ A, π b ∈ B and so forth. By distinct, we mean for example that only one of the equivalent processes π 1 π 2 → π 3 π 4 and π 2 π 1 → π 3 π 4 appears in P . It is then possible to define an effective amplitude square as: This amplitude can then be used to compute an effective cross section σ P (s) using the standard procedure. Of course, factors for identical incoming or outgoing particles should not be reintroduced as they are already taken into account. No average on incoming particles is taken either. One can then compute the thermally-averaged cross section

A.3 3 → 2 scattering
Consider a 3 → 2 scattering process p given by p : π a π b π c → π d π e . (A. 19) In the limit that all pions are degenerate in mass, the corresponding cross section is [4] σv 2 p = where sign({a, b, c, d, e}) is the sign of the permutation of {a, b, c, d, e}, e.g the sign of {1, 2, 3} is +1 and that of {1, 3, 2} is −1. Note that eq. (A.4) implies that five different pions must be involved in the Lagrangian term. This means that there is no subtlety concerning identical particles in the ingoing or outgoing states. A complex pion can however appear as both itself and its conjugate, meaning that a pion can appear on one side of a 3 → 2 process and the other at the same time.
We neglect the mass splitting between the pions for two reasons. First, the 3 → 2 processes that we consider are always between particles that are close in mass and have far more mass incoming than outgoing. The mass splitting is therefore a small effect, specially considering that these processes take place at high temperature. Second, in the end, the dark matter abundance is mainly determined by when the stable pions decouple from the unstable ones. The 3 → 2 processes do not affect this much.
As was done for 2 → 2 processes, an effective cross section can be defined by merging processes involving pions from the same multiplets. Define P : π A π B π C → π D π E as the set of distinct p : π a π b π c → π d π e such that π a ∈ A, π b ∈ B and so forth. The effective cross section is

A.4 Boltzmann equation
The results of the previous sections can then be combined to write down Boltzmann equation. Let's begin by introducing additional notation: • Define P 2→2 as the set of all distinct 2 → 2 scatterings between representations.
• Define P 3→2 as the set of all distinct 3 → 2 scatterings between representations.
• Define N A as the size of the representation A.
• Define the entropy density as s.
• Define H as the Hubble constant.
• DefineŶ A as the number density of pions from A per entropy density divided by N A .
• DefineŶ eq A as the equilibrium value ofŶ A .
• Define #A P in as the number of pions from A in the in state of a process P .
• Define #A P out as the number of pions from A in the out state of a process P .
• Define Γ A as Γ A 1/γ A , where 1/γ A is the average of the inverse of the gamma factor for a pion from multiplet A.
With all this notation established, the Boltzmann equation can be written as: σv P :π P 1 π P 2 →π P 3 π P 4 + s eq P 5Ŷ P 4Ŷ P 5 σv 2 P :π P 1 π P 2 π P 3 →π P 4 π P 5 where the first sum is over processes such that #A P in > #A P out to avoid double counting.

A.5 Numerical procedure
We now proceed to describe how we solve the Boltzmann equation in practice. At first, the pions are kept at their thermal equilibrium value by 3 → 2 processes. The validity of this assumption is tested periodically by comparing the rate of 3 → 2 processes to the Hubble constant. When this approximation becomes remotely close to failing, the pion densities are evolved as a whole by assuming that they are in chemical equilibrium with each other. The validity of this assumption is tested periodically by computing the rates of 2 → 2 scattering if one of the pion densities were to vary a little and comparing this rate to the Hubble constant. When one particle comes remotely close to failing this assumption, we begin to evolve it individually. The densities of the other pions are computed as before but also by keeping track of how many were converted to pions that are not in chemical equilibrium anymore. Once none of the pions are in chemical equilibrium anymore, all pions are evolved using the full Boltzmann equation.

B Symmetry breaking benchmarks
In this appendix, we present relevant information for each benchmark scenario of symmetry breaking of table 1.

B.1 Notation
First, we introduce some notation for the pions. The pion matrix can be written in all cases we consider as

JHEP02(2020)196
The matrix C is Hermitian, traceless and can be written for all patterns of chiral symmetry breaking with N flavors as: The matrix D can be written for each pattern of chiral symmetry breaking as:

(B.3)
The generalization to larger N is trivial in both cases. The different types of pions are: • Pions π d i consist of a linear superposition of a quark and the corresponding antiquark. They are present for all patterns a chiral symmetry breaking.
• Pions π o i consist of a quark and a different antiquark. They are present for all patterns a chiral symmetry breaking.
• Pions π a i consist of two identical quarks. They are present for SU(2N ) → SO(2N ) only.
Obviously, pionsπ j i are the conjugate of π j i . The generators t a can then be read by simply decomposing Π = t a π a . They are normalized such that Tr[t a t b ] = δ ab , i.e. 1 if π a is the conjugate of π b and 0 otherwise.

B.2 Benchmarks
We then proceed to discuss the benchmarks. For each benchmark, we will provide the list of representations of h, the pions that constitute them, their masses and their stability. They are ordered from heaviest to lightest. In some cases, the mass eigenstates are linear -32 -

JHEP02(2020)196
combinations of several pions of eq. (B.1). We express such states as the combination of these pions within brackets, eg. (π d 1 , π d 2 ). A pion will be considered stable if there does not exist a set of pions which can be combined to reproduce the same charges and whose total mass is inferior. We will also indicate in the column ID whether that representation contributes to indirect detection or not. If it does, we provide an example. Assumed relations between the quark masses are included. Relevant comments are provided. When multiple U(1) symmetries are present, they can be combined in different ways and we provide one example only. Their charges are normalized to the smallest integers. For representations of h 1 × h 2 with h 1 non-Abelian and h 2 Abelian, the notation (n, ±m) stand for the combination of (n, +m) and (n, −m). Generalizations of this notation are straightforward. • m 1 > m 2 > m 3

Comment(s):
• Any other hierarchy of dark quark masses would lead to this symmetry breaking pattern being part of category I as long as no two dark quark masses are equal.
• Two relations are relevant: • Any other hierarchy of dark quark masses would lead to this symmetry breaking pattern being part of category II as long as no two dark quark masses are equal.