Cosmological Bounds on Non-Abelian Dark Forces

Non-Abelian dark gauge forces that do not couple directly to ordinary matter may be realized in nature. The minimal form of such a dark force is a pure Yang-Mills theory. If the dark sector is reheated in the early universe, it will be realized as a set of dark gluons at high temperatures and as a collection of dark glueballs at lower temperatures, with a cosmological phase transition from one form to the other. Despite being dark, the gauge fields of the new force can connect indirectly to the Standard Model through non-renormalizable operators. These operators will transfer energy between the dark and visible sectors, and they allow some or all of the dark glueballs to decay. In this work we investigate the cosmological evolution and decays of dark glueballs in the presence of connector operators to the Standard Model. Dark glueball decays can modify cosmological and astrophysical observables, and we use these considerations to put very strong limits on the existence of pure non-Abelian dark forces. On the other hand, if one or more of the dark glueballs are stable, we find that they can potentially make up the dark matter of the universe.


I. INTRODUCTION
New gauge forces may be realized in nature beyond the SU (3) c × SU (2) L × U (1) Y structure of the Standard Model (SM). If a new gauge force connects directly to SM matter, it must have a characteristic mass scale above about a TeV to be consistent with experimental tests of the SM [1][2][3]. On the other hand, new dark gauge forces that couple only very weakly to the SM can be significantly lighter [4][5][6]. Such dark forces can be very challenging to probe directly in experiments, and in many scenarios the strongest bounds on them come from astrophysical and cosmological observations [7][8][9][10].
In this work we investigate the cosmological evolution and constraints on new non-Abelian dark forces. Such dark forces are well-motivated, and arise in string theory constructions [11], models of dark matter or baryogenesis , neutral naturalness scenarios [33][34][35][36][37], and within the hidden valley paradigm [38][39][40]. The requirement of gauge invariance in theories of non-Abelian dark forces implies that the new gauge vector bosons can only couple to the SM through non-renormalizable operators [39,40]. This stands in contrast to Abelian dark forces that can connect to the SM at the renormalizable level through kinetic mixing with hypercharge. As a result, direct low-energy searches for non-Abelian dark forces are very difficult, and cosmological observations usually provide the most powerful tests of them [19][20][21][22][23][24][25][26][27][28][29][30][31][32].
The particle spectrum in theories of non-Abelian forces is diverse and complicated, and depends on both the gauge group and the representations of the matter fields charged under it. We focus on the minimal realization of a non-Abelian dark force consisting of a pure Yang-Mills theory with a simple gauge group G x . Such theories can be described in terms of self-interacting dark gluons at high energies, but are expected to confine to a set of gauge-neutral glueball bound states below a dynamically generated confinement scale Λ x [41]. Both phases can be realized in the hot early universe, with a transition from the gluon phase to the glueball phase occurring as the temperature (of the dark sector) falls below the critical temperature T c ∼ Λ x [42].
If the visible and dark sectors do not interact, they evolve independently with distinct temperatures T and T x . After confinement at T x = T c , the dark glueballs undergo a complicated freezeout process. The energy density of the dark sector is dominated by the lightest glueball state, which on general grounds is expected to have J P C = 0 ++ [43]. The lightest 0 ++ number density changes mainly through (3 → 2) self-annihilation processes [25,44]. While these reactions are active, the dark temperature changes very slowly, only falling off as the logarithm of the cosmological scale factor [45,46]. As a result, the lightest glueballs form a massive thermal bath in which the other heavier glueballs annihilate through 2 → 2 processes and eventually freeze out [44,47,48]. In the end, a collection of relic glueball densities is left over, dominated by the 0 ++ with exponentially smaller yields for the heavier states [44,48].
The process of glueball freezeout can change drastically if there are operators that connect the visible and dark sectors. Such operators are always expected at some level; quantum gravitational effects are thought to induce gauge-invariant operators involving both SM and dark sector fields suppressed by powers of the Planck mass [49][50][51][52]. Even stronger connections can arise if there exist new matter fields that couple directly to both the visible and dark sectors [39,40]. As long as the new physics generating these operators is much larger than the confinement scale, their effects can be parametrized in terms of a set of non-renormalizable connector operators.
With connectors, energy can be transferred between the dark and visible sectors [14,25,29,30,44]. After confinement, connector operators can also modify the glueball freeze-out dynamics and induce decays of some or all of the dark glueballs to the SM. If one of the glueballs is longlived or stable, it will contribute to the density of dark matter (DM) [21]. However, glueball lifetimes that are not exceedingly long will inject energy into the cosmological plasma and modify the standard predictions for big bang nucleosynthesis (BBN) [53,54] and the cosmic microwave background (CMB) [55,56], as well as act as astrophysical sources of cosmic and gamma rays [9].
The aim of this work is to estimate the bounds on pure non-Abelian dark forces in the presence of connector operators from cosmology and astrophysics. We focus mainly on the dark gauge group G x = SU (3) with glueball masses above m 0 ≥ 100 MeV, and we study the leading connector operators between the dark vector bosons and the SM with characteristic mass scale M m 0 . As an initial condition, we assume inflation (or something like it) followed by preferential reheating to the visible sector to a temperature above the confinement scale but below that of the connectors. With these assumptions, we find very strong limits on non-Abelian dark forces.
Cosmological effects of dark gluons and glueballs were studied previously in Refs. [9,14,19,21,22,24,25,[29][30][31]44]. We extend these earlier works with a more detailed analysis of the leading (2-body) connector operators and their effects on energy transfer between the visible and dark sectors. We also investigate the effects of heavier glueballs in the spectrum beyond the lightest mode, and we show that the lightest C-odd glueball can play an important role in some cases and even make up the observed DM density when it is long lived or stable.
Following this introduction, we discuss the general properties of glueballs relevant to our analysis in Sec. II. Next, we present the leading connector operators to the SM and investigate their implications for glueball decays in Sec. III. In Sec. IV we study the cosmological evolution of the dark gauge theory and we compute glueball yields both with and without connector operators. These results are then applied to derive cosmological constraints on dark glueballs in Sec. V. Finally, Sec. VI is reserved for our conclusions. Some technical details about gluon thermalization and the cosmological and astrophysical bounds we apply are collected in Appendices A and B.

II. GLUEBALL PROPERTIES
Glueballs have been studied using a variety of methods for a wide range of non-Abelian gauge groups [41,57]. In this section we review and derive some general results for SU (N ) glueballs that will be essential for the analysis to follow.

A. Glueball Masses
Detailed lattice studies of glueballs have been performed for SU (N ) with N = 2, 3, . . ., and a number of stable states are found. Since the minimal Yang-Mills action respects angular momentum (J), parity (P ), and charge conservation (C), glueballs can be classified according to their J P C quantum numbers. The lightest state is found to have J P C = 0 ++ [58,59], as expected on general grounds [43]. The masses and quantum numbers of the stable glueballs found for SU (2) and SU (3) are listed in Tab. I. They are expressed in terms of the length scale r 0 where the inter-gluon potential goes from Coulombic to linear, and for SU (3) is related to the strong coupling scale by r 0 Λ x = 0.614(2)(5) [60].
For reasons to be explained below, we focus our attention on two specific glueball states (for SU (N ≥ 3)): the lightest overall J P C = 0 ++ glueball in the spectrum, together with the lightest C-odd glueball with J P C = 1 +− . With the gauge group SU (3), the mass of the lightest 0 ++ glueball is m 0 6.9 Λ x , and the 1 +− mass is m 1 = 1.71(5) m 0 [59].
Going beyond SU (3) to larger N , the glueball mass spectrum is found to be similar, with mass corrections suppressed by powers of 1/N 2 [61]. Similar results are also expected for other gauge groups with non-vanishing anomaly coefficient d abc = tr(t a {t b , t c }), where t a is the generator of the fundamental representation [41] (which we normalize according to tr(t a t b ) = δ ab /2 for the N of SU (N )). However, let us point out that for SU (2) and other groups with vanishing d abc such as SO(2N + 1) and Sp(2N ), there are no C-odd states in the spectrum [39]. A further extension J P C m r 0 (N = 2) m r 0 (N = 3)  of the minimal Yang-Mills theory is the inclusion of a topological theta term. This would break P and T , but not C. It would also shift the glueball masses [62], and induce mixing between glueball states with different P quantum numbers [62,63].

B. Glueball Couplings and Matrix Elements
Glueball self-couplings and transition matrix elements are needed to compute their cosmological evolution. These quantities have not been studied in as much detail on the lattice as the glueball mass spectrum. Here, we collect the relevant existing lattice results, and we use naive dimensional analysis (NDA) [64][65][66] and large-N scaling [67,68] to make estimates when no lattice data is available.
Glueball interactions are expected to be perturbative in the limit of large N (for an underlying SU (N ) gauge group), and this motivates writing an effective Lagrangian in terms of glueball fields. Combining the N scaling of gluon n-point functions with dimensional analysis suggests the form where φ represents a glueball field interpolated by a single-trace gluon operator, m x is a characteristic glueball mass scale, and F (x, y) is a smooth function that is finite as N → ∞. Expanding this function in a power series and rescaling to obtain a canonical kinetic operator, the effective Lagrangian becomes where the coefficients a n are expected to be of order unity. Note that shifting the gluon field to remove the linear term does not alter this general form. In the analysis to follow, we identify m x = m 0 with the mass of the lightest glueball. We will also need glueball matrix elements in the analysis to follow. Specific glueball states can be identified with gauge invariant gluon operators, in the sense that the operators can create one-particle glueball states from the vacuum. For example [39], Here, X µν = X a µν t a is the dark gluon field strength contracted with the generators of the fundamental representation of the group normalized to tr(t a t b ) = δ ab /2.
The two matrix elements of greatest interest to us are where the estimates on the right hand sides are based on large-N and NDA, and α x = g 2 x /4π is the dark gauge coupling. In the second line, we have also suppressed the Lorentz structure of the matrix element, µναβ p α ε β , where p α is the outgoing momentum and ε β is the polarization of the initial state [39]. The first of these matrix elements, F S 0 ++ , has been computed on the lattice for N = 3 with the result [59,69] which agrees reasonably well with our large-N and NDA estimate and is scale independent. In contrast, the second matrix element has not been calculated on the lattice. We use the lattice value of F S 0 ++ and the NDA estimate α 3/2 x M 1 +− 0 ++ = 4π/N m 3 x in the analysis to follow.

III. CONNECTIONS TO THE SM AND GLUEBALL DECAYS
With the SM uncharged under the dark gauge group G x , gauge invariance forbids a direct renormalizable connection of the dark gluons to the SM. However, massive mediator states that couple to both sectors can generate non-renormalizable operators connecting them. If the characteristic mass scale of the mediators is M Λ x , the leading operators have mass dimension of eight and six, and take the form [39,40] where X and F SM refer to the dark gluon and SM field strengths. If present, these operators allow some or all of the glueballs to decay to the SM. In this section we illustrate simple mediator scenarios that generate these operators, and we compute the glueball decay rates they induce.

A. Dimension-8 Operators
Dimension-8 operators of the form of Eqs. (7,8) lead to glueball decays with characteristic rate Here, we present an explicit scenario of mediator fermions that generates these operators and we compute the glueball decay rates they induce. Before proceeding, it is helpful to organize the dimension-8 operators according to a dark charge conjugation operation C x under which X a µ → −η(a)X a µ , where η(a) is the sign change of the fundamental generator t a under charge conjugation [70], with the SM vector bosons being invariant. The operators of Eq. (7) are even under C x and those of Eq. (8) are odd. Furthermore, C x coincides with the C x -number assignments of the glueball states. Correspondingly, the operators of Eq. (7) only allow direct decays of C x -even glueballs to the SM, or glueball transitions from even to even or odd to odd. In particular, at d = 8 the operator of Eq. (8) is required for the lightest C x -odd 1 +− glueball to decay.
Consider now a set of massive vector-like fermions with masses M r ∼ M Λ x , each transforming as a fundamental or antifundamental under G x = SU (N ) and the representation r of the SM gauge group (defined with respect to the left-handed component of the fermion). Direct collider and precision electroweak limits on such fermions imply M r 100 GeV if they only have electroweak charges, and M r 1000 GeV if they are charged under QCD [39,40,71]. The effective Lagrangian generated by integrating the fermions out is [39] Here, the dark gluon operators S, P, and Ω (1,2) µν correspond to Eq. (3), and the coefficients χ i are given by where the sums run over the SM representations r of the fermions, and ρ r = M r /M . For each such representation, we define sub-representations r = (r 1 , r 2 , r 3 ) with respect to the SM gauge The quantity d(r i ) is the number of copies of the i-th subrepresentation within r, and T 2 (r i ) is the trace invariant for that factor (normalized to 1/2 for the N of SU (N ) and Y 2 for U (1) Y ). 1 Generic representations of mediator fermions break the dark charge conjugation number C x explicitly and generate both operator types of Eqs. (7,8). This is explicit in Eq. (11), with both even (χ i = 0) and odd operators (χ Y = 0). However, there exist mediator fermion combinations that preserve C x [72] and yield χ Y = 0. From Eq. (14), we see that this requires a specific combination of fermion charges as well as masses. The presence of masses also implies that C x can be broken softly. In contrast, the χ i coefficients of Eq. (13) are positive semi-definite and not subject to cancellation.
The C x -preserving operator of Eq. (11) allows direct decays of the 0 ++ glueball to pairs of SM vector bosons. The corresponding decay widths are [39] where m 0 = m x is the 0 ++ glueball mass, F S 0 ++ is given by Eq. (4), N 2 c − 1 = 8, the χ i are defined in Eq. (13), with s W being the sine of the weak mixing angle. Note that the decay width to gluons in Eq. (15) only applies for m 0 1 GeV; at lower masses the final states consist of hadrons. We do not attempt to model this hadronization, and instead we apply a factor of 1 − (2m π /m 0 ) 2 to the decay width. In evaluating the width of Eq. (15), we take α 3 at scale m 0 since the corresponding gluon operator is renormalized (at one-loop) in the same way as the standard field strength operator.
Decays of the lightest 1 +− glueball occur through the C x -odd operator term in Eq. (11), with the leading decay channels expected to be 1 +− → 0 ++ + {γ, Z}. The widths are [39] with m 1 = m 1 +− , and M 1 +− 0 ++ defined in Eq. (5). The total decay lifetimes τ = 1/Γ of the 0 ++ and 1 +− glueball states from the dimension-8 operators above with χ i = χ Y = 1 and G x = SU (3) are shown in the left and right panels of Fig. 1. In the upper left of both plots, we mask out the regions with m 0 > M/10 where our treatment in terms of effective operators breaks down. The dotted, solid, and dashed lines indicate reference lifetimes of τ = 1/Γ = 0.1 s, 5 × 10 17 s, 10 26 s. These lifetimes correspond to decays that occur early in the history of the Universe, at the present day, and long lived glueballs, respectively. Both decay rates follow the approximate scaling of Eq. (10). All other known (SU (3)) glueballs can decay through these dimension-8 operators as well with parametrically similar rates, although there can be numerically significant differences due to coupling factors and phase space [39].

B. Dimension-6 Operators
Glueball decays through the dimension-6 operator of Eq. (9) proceed with characteristic rate We present here two mediator scenarios that generate the operator of Eq. (9) and we compute the decay rates they induce. Our first mediator scenario follows Ref. [40] and consists of mediator fermions with Yukawa couplings to the SM Higgs boson. A minimal realization contains a vector-like SU (2) L doublet P with gauge quantum numbers (r x , 1, 2, −1/2), and a vector-like singlet N with quantum numbers (r x , 1, 1, 0) together with the interactions [40,71] For M N , M P m h , the leading glueball effective operator from integrating out the fermions can be obtained using the low-energy Higgs theorem [73], where M 2 M P M N and T 2 (r x ) = 1/2 is the trace invariant of the fermion representation r x under the dark gauge group G x . In addition to the dimension-6 operator above, the massive fermions also generate dimension-8 operators of the form of Eq. (11). A second mediator scenario consists of a complex scalar Φ x charged under the dark gauge group with a Higgs-portal coupling, Applying the low-energy Higgs theorem to this state (for M Φ m h ), we find In passing, we note that the Higgs portal coupling of Eq. (25) respects dark C x number. The operator generated in either mediator scenario can be written in the form with the dimensionless coefficient y ef f . Since this operator is even under C x , it only allows direct decays of C x -even glueballs to the SM, or even-to-even or odd-to-odd glueball transitions. It was shown in Ref. [40] that this is sufficient to allow all known SU (3) glueballs to decay, except for the 1 +− and 0 −+ modes. The absence of a 1 +− decay follows from C x considerations, while the conclusion for 0 −+ is a result of spin and parity, rather than C x . This mode can decay at the dimension-6 level if a topological dark gluon term is added to the UV Lagrangian or by extending to a two-Higgs doublet model [40].
Using the parametrization of Eq. (27), the direct decay of the 0 ++ glueball to the SM has rate [40] where √ 2 H = 246 GeV is the electroweak vacuum expectation value, F S 0 ++ is defined in Eq. (4), m h = 125 GeV is the Higgs mass, Γ h = 4.1 MeV is the Higgs width, and Γ h (m h → m 0 ) is the total width the SM Higgs would have if its mass were m 0 (and includes decays to Higgs final states for m 0 > 2m h ). We evaluate this width using the expressions of Refs. [74,75].
In Fig. 2 we show the decay lifetime τ = 1/Γ of the 0 ++ glueball from the dimension-6 (and dimension-8) operators above with y ef f = 1 and G x = SU  Fig. 1, we see that it is parametrically long-lived compared to the 0 ++ when both dimension-6 and dimension-8 operators are present.

C. Decay Scenarios
Based on the discussion above, we present four glueball decay scenarios organized by the dimensions of the relevant decay operators and the dark conjugation charge C x : 1. Dimension-8 decays with broken C x In this scenario glueballs decay exclusively through the dimension-8 operators of the form of Eq. (11). All glueballs are able to decay with parametrically similar rates. To realize this scenario, we use the effective interactions in Eq. (11) with χ i = χ Y = 1.

Dimension-8 decays with exact C x
This scenario is similar to the first, but now with χ Y = 0. Conservation of C x implies that the lightest 1 +− glueball is stable. The other glueballs are all able to decay with parametrically similar rates.
3. Dimension-6 decays with broken C x Glueball decays occur through the dimension-6 operator of Eq. (27) and the dimension-8 operators of Eq. (11). We realize the scenario by setting y ef f = 1 together with χ i = χ Y = 1. With the exception of the 1 +− mode (and possibly the 0 −+ ), glueballs decay primarily through the dimension-6 operator. In contrast, the 1 +− glueball only decays through the C x -breaking dimension-8 operator with a parametrically suppressed rate, making it much longer-lived than the other glueballs, which in turn leads to different cosmological scenarios when considering the constraints we can place on this model.

Dimension-6 decays with exact C x
Decays occur through the dimension-6 operator of Eq. (27) and the C x -conserving terms in Eq. (11). We realize the scenario by taking y ef f = 1, χ i = 1, and χ Y = 0. The 1 +− glueball is stable, while the other glueballs decay mainly through the dimension-6 operator.
We study the cosmological implications of these four decay scenarios in the analysis to follow.

IV. GLUEBALL DENSITIES IN THE EARLY UNIVERSE
Glueballs are formed in the early universe in a confining transition as the dark sector temperature T x falls below a critical temperature T c ∼ m 0 . After they are created, the glueballs undergo a complicated freezeout process involving a range of 2 → 2 and 3 → 2 reactions. These dynamics become even more complicated when the dark sector connects to the SM through the operators discussed above, with new effects such as energy transfer between the visible and dark sectors and glueball decays. In this section we review the formation and freezeout of glueballs in the absence of connectors to the SM, and we investigate how this picture changes when connectors are present.

A. Glueball Formation and Freezeout without Connectors
In the absence of operators that connect to the SM, the visible and dark sectors do not thermalize with each other. We assume that enough energy is liberated by reheating following primordial inflation (or something similar) that both sectors are able to thermalize independently with temperatures T and T x [76], and furthermore that T x ≥ T c at this point. 2 As the universe expands and cools, dark glueballs are formed in a confining transition. This transition has been studied in detail using lattice methods for G x = SU (N ) [42,[77][78][79][80], and the critical temperature of the transition is found to be [79] T c r 0 = 0.709(6) + 0.546 (22) corresponding to T c m 0 /5.5 for N = 3. 3 The confining transition is also found to be secondorder for N = 2, very weakly first-order for N = 3, and increasingly strongly first-order for larger N [42,78]. Generalizing the analysis of Ref. [24] as in Ref. [44], we estimate that the fractional entropy change in the transition is negligible for N 10 over the full range of parameters considered in this work. Following the transition, as T x falls below the critical temperature, the glueball masses quickly settle to their zero-temperature values [82,83]. Based on these results, we assume the phase transition occurs instantaneously at T x = T c . This should be a good approximation for N = 3, but could be inaccurate for much larger N .
Entropy is conserved independently in both sectors while kinetic equilibrium is maintained. This implies that the ratio of entropy densities s and s x in the two sectors remains constant, We take R as an input to our calculation; in the absence of connectors its value is set by the unspecified dynamics of reheating after inflation [76]. However, we do assume R < 1 corresponding to preferential reheating to the visible sector. For T x T c and G x = SU (N ), the entropy ratio is related to the temperatures in the two sectors by where g * S is an effective number of degrees of freedom in the visible sector at temperature T . This ratio will be maintained through the confining transition provided it is not too strongly first order [24]. Once formed, dark glueballs interact with each other and undergo a freezeout process in which they depart from thermodynamic equilibrium and develop stable relic densities. This process was studied in detail in Refs. [24,25,44]. In the last work, the evolution of glueball numbers was computed numerically using a network of Boltzmann equations containing the most important 2 → 2 and 3 → 2 reactions, with thermally averaged cross sections estimated using the glueball effective Lagrangian of Eq. (2).
For the purposes of our cosmological analysis of glueball effects to follow, the results of Ref. [44] are captured to an excellent approximation by a simplified two-state model for the densities of the 0 ++ and 1 +− glueballs. In this model, the evolution equations for the 0 ++ density n 0 and the 1 +− density n 1 are [44] dn 0 dt where H is the Hubble factor,n i are the equilibrium number densities at temperature T x , and σ 32 v and σ 22 v correspond to the 3(0 ++ ) → 2(0 ++ ) and (1 +− 1 +− ) → (0 ++ 0 ++ ) processes. In detail, the Hubble factor is given by where ρ is the energy density of the SM and ρ x is that of the dark sector. Since kinetic equilibrium is expected to hold in the dark sector throughout the freezeout process, the number densities take the form where E i = m 2 i + p 2 , g i is the number of internal degrees of freedom, µ i is a chemical potential, and T x is the common dark sector temperature. The precise definition of the equilibrium densities is thenn i = n i (T x , µ i → 0). For the thermally averaged cross sections, we estimate them using Eq. (2): To evaluate Eqs. (32,33), it is necessary to track the time evolution of the dark temperature T x . This can can achieved using the constancy of the entropy ratio R together with Prior to freezeout of the 0 ++ mode (and after dark confinement), the dark temperature falls as T x ∝ 1/ ln(a) due to the energy injected by 3 → 2 annihilations [45]. After 0 ++ freezeout, the dark temperature falls as T x ∝ a −2 .
Our two-state simplified model provides an excellent approximation of the full analysis of Ref. [44]. A central feature of the analysis is that near equilibrium at T x < T c the 2 → 2 reactions are parametrically faster than the 3 → 2. This keeps the ratio of the 0 ++ and 1 +− densities close to the equilibrium ratio throughout the freezeout process, with ∆m = (m 1 − m 0 ) and x x = m 0 /T x . Thus, the 2 → 2 reactions push the 1 +− density to be exponentially smaller than the 0 ++ . The exponential suppression of heavier modes also means that the freezeout of the 0 ++ mode can be computed reliably in isolation, neglecting the effects of the 1 +− and keeping only the 3 → 2 reactions. Calculations of single-state freezeout through 3 → 2 annihilation have been performed in Refs. [25,45,47,48]. For freezeout at Numerically, we find x f o x ∈ [5,20] for R ∈ [10 −12 , 0.1] and m 0 ∈ [10 −3 , 10 9 ] GeV. The less-abundant 1 +− mode freezes out in the background of the massive bath of 0 ++ glueballs [44,48]. This occurs after 0 ++ freezeout, but before the kinetic self-equilibration of the 0 ++ states is lost.
In Fig. 3  densities are set by the non-perturbative dynamics of the confining phase transition. As expected, the 1 +− yield is always much lower than the 0 ++ yield.
Going beyond the two-state model, our arguments regarding the exponential suppression of the 1 +− density relative to the 0 ++ also apply to the other heavier glueball modes [44]. The total glueball relic density is strongly dominated by the 0 ++ density, while 2 → 2 annihilation reactions push the heavier glueball densities to much smaller values. In fact, these reactions tend to be much more efficient for the other heavier glueballs than the 1 +− due to coannihilation with the 0 ++ . For example, 2 ++ + 0 ++ → 0 ++ + 0 ++ efficiently depletes the second-lightest 2 ++ glueball up to very large x x , while reactions such as 1 −− + 0 ++ → 1 +− + 0 ++ quickly transfer the density of heavier C-odd glueballs to the lighter 1 +− . Conservation of C x number in the dark sector implies that the net density of C-odd glueballs cannot be reduced by coannihilation. As a result, the 1 +− state generally develops the second largest relic density, with the densities of the other dark glueballs being much smaller. This, combined with the unique decay properties of the 1 +− glueball when connectors are included, is the reason why we only consider the effects of the 0 ++ and 1 +− glueballs in our analysis of glueball cosmology.

B. Glueball Freezeout with Connectors
Connector operators can modify the freezeout of glueballs in a number of ways. Scattering and decay reactions mediated by such operators transfer energy between the visible and dark sectors, and may allow them to thermalize. Decays through the connector operators after confinement also deplete glueballs, and can occur before or after the freezeout of the various (3 → 2) and (2 → 2) reactions. We investigate these effects here, both before and after confinement, with a focus on the 0 ++ and 1 +− glueballs. Our goal is to compute the yields of these species prior to their decay.
As in the freezeout analysis without connectors, we take as an initial condition primordial inflation (or something like it) with preferential reheating to the visible sector characterized by a temperature T RH that is larger than the confinement transition temperature T c m 0 /5.5. With connectors, we also assume T RH M . Reheating above the connector scale M is likely to thermalize the dark and visible sectors at T RH , and can produce a relic abundance of the connector particles themselves. These can have interesting cosmological effects in their own right, acting as quirks if they carry G x charge [34,85,86], and potentially creating dark glueballs nonthermally [31,32,87,88]. By taking T RH M , the production of connector particles in the early universe is strongly suppressed allowing us to focus on the effects of the glueballs.

Energy Transfer before Confinement
Consider first the transfer of energy at temperatures T well above the confinement temperature T c . In the absence of connectors, preferential reheating to the visible sector produces T x T . Connector operators allow reactions of the form SM + SM ↔ X + X that transfer energy from the visible sector to the dark sector. For T x > T c , the evolution equation for the energy density of the dark sector is [89,90] where ∆E · σv is the thermally averaged energy transfer cross section for X + X → SM + SM, n x is the dark gluon number density, and n x =g x (ζ(3)/π 2 )T 3 is the value it would have in full equilibrium with the visible sector withg x dark gluon degrees of freedom (equal tog x = 2(N 2 −1) for G x = SU (N )). 4 For T x T , the n 2 x term on the right side above dominates and leads to a net energy transfer to the dark sector. This transfer saturates and ceases when T x → T and n x → n x .
For visible radiation domination with constant g * , Eq. (41) can be rewritten as d dT With the connector operators of Eqs. (7,9) and T T x , the right side of Eq. (42) takes the parametric form where n = 4, 8. Integrating from temperature T to the reheating temperature T RH , the approximate solution is This expression is dominated by the contribution near the reheating temperature, and represents the contribution to the dark energy density from transfer reactions. The approximate forms of Eqs. (44,45) are only valid for T < T RH and T > T x ≥ T c . The first of these conditions corresponds to the upper limit on the era of radiation domination. An even higher radiation temperature can be achieved prior to reheating, but for standard perturbative reheating and n < 29/3 9.67 we find that the energy transfer before the radiation era is also dominated by reactions near T ∼ T RH . The second condition T > T x ≥ T c is needed to justify our neglect of the n 2 x term on the right side of Eq. (42) and our assumption of a deconfined phase. As T x approaches T due to the energy transfer, this term becomes important and the net energy transfer goes to zero, corresponding to the thermalization of the two sectors.

Motivated by these considerations, let us define
This represents the contribution to the dark sector energy from thermal transfer in the vicinity of temperature T . Thermalization occurs when whereg x is the number of dark gluon degrees of freedom. Let T th be the temperature that solves Eq. (48) as an equality. If T th < T c , the visible and dark sectors remain thermalized at least until confinement. Conversely, if T th > T c thermalization is lost at T = T th and the dark and visible sectors evolve independently thereafter with separately conserved entropies. The dark to visible entropy ratio R is constant for T < T th and depends on reheating. If T th < T RH , thermalization occurs after reheating and is maintained until T = T th . The entropy ratio R (for T th > T c ) after thermalization ceases is then Thermalization need never have occurred after reheating if T RH < T th . In this case, (for T th > T c ) we can define This implies a lower bound on the entropy ratio of In general, lower reheating temperatures allow for smaller values of R. We define R min to be the value of R such that T x RH = T c , the lowest possible reheating temperature given our assumption of T x RH ≥ T c . 5 When T th > T c , the range of R values is therefore R min ≤ R ≤ R max . In Appendix A we present explicit expressions for the collision term ∆C needed to compute the energy transfer ∆(ρ x /T 4 ) via Eq. (46). The results obtained for R min are shown in Fig. 4 in the In the cyan region in the right panel, thermalization between the visible and dark sectors is maintained at least until confinement.

Evolution of the 0 ++ Density
Glueballs form at T x = T c and undergo freezeout, transfer, and decay reactions. In the absence of connectors, the dominant glueball species is the lightest 0 ++ mode. To track its evolution with connector operators, it is convenient to organize the analysis according to the thermalization temperature T th , computed above in the unconfined phase, relative to the confinement temperature.
T th < T c : This condition implies that thermalization is maintained at least until confinement, and thus we expect T = T x = T c as an initial condition for the glueball evolution. To compute the 0 ++ density and thermal transfer after confinement we adapt the analysis of Refs. [93,94] based on Refs. [89,90,95], which is applicable here since T, T x ≤ T c m 0 /5.5. If thermal equilibrium is maintained independently within both the dark and visible sectors, the dark temperature evolves as [93,94] where C ρ and C n are the collision terms appearing in the evolution equations for the 0 ++ energy and number densities. The Hubble term in Eq. (52) gives the usual 1/a 2 redshifting of the effective temperature of an independent massive species, while the second term describes energy transfer from scattering and decay processes. The explicit forms of the collision terms are wheren 0 = n 0 (T x ) and n 0 = n 0 (T ), as well as The only new piece in these expressions is the elastic scattering term n 0 n SM σ el v·∆E in Eq. (54). It corresponds to reactions of the form SM + 0 ++ → SM + 0 ++ , and was studied in detail in Refs. [89,90].
Combined in Eq. (52), the (3 → 2) scattering term from Eq. (53) tends to heat the dark glueballs, and the elastic scattering and decay terms tend to drive T x → T . Applied to the 0 ++ glueball with either the dimension-8 or dimension-6 connector operators, we find that thermalization below confinement implies Γ 0 > H(T = m 0 ). Thus, the 0 ++ density simply tracks the equilibrium value with temperature T following confinement. 6 T th > T c : With T th > T c , the visible and dark sectors are not thermally connected at confinement, and thus T ≥ T x at this point with a well-defined entropy ratio in the range R min ≤ R ≤ R max . Using the scaling arguments applied above, it can be shown that R ≥ R min implies T ≤ m 0 when the 0 ++ decays set in at Γ 0 H(T ). 7 The evolution equations for the 0 ++ number density and temperature can thus be written as (to leading order in T x /m 0 ) where againn 0 is the equilibrium value at temperature T x and n 0 is the equilibrium value at temperature T . Note that we have neglected the elastic scattering term because it can be shown to be parametrically small relative to the Hubble term for T < T th and R ≥ R min . When the decay terms are neglected, the evolution equations of Eqs. (55,56) are equivalent to those we used previously with no connector operators (to leading order in T x /m 0 ). Glueball decays only become significant when Γ 0 H(T ), and quickly drive T x → T and n 0 → n 0 . It follows that our previous analysis without connectors can be applied to compute the 0 ++ relic yield prior to decay (which may occur before freezeout). The only significant effect of energy transfer on this calculation is to limit the range of the initial entropy ratio to R min ≤ R ≤ R max .

Evolution of the 1 +− Density
Even though the lightest 0 ++ glueball dominates the total glueball density and controls the dark temperature prior (and even after) its decay, the heavier 1 +− glueball can also be relevant for cosmology due to its longer lifetime. Recall that the 1 +− is parametrically long-lived relative the 0 ++ in the decay scenarios 2-4 listed in Sec. III C, where the 0 ++ decays through a dimension-6 operator while the 1 +− is stable or only decays at dimension-8. Even in decay scenario 1, where both states decay at dimension-8, the 0 ++ decay rate tends to be larger than the 1 +− by a factor of (N 2 c − 1)α 3 /α. The evolution of the 1 +− density is sensitive to the 0 ++ density in several ways. Prior to decay, the 0 ++ density acts as a massive thermal bath that cools very slowly relative to the visible temperature, thereby delaying the freezeout of the 1 +− state. This thermal bath collapses and disappears when the 0 ++ decays, which can hasten 1 +− freezeout. If the 0 ++ density is large when it decays, the entropy transferred to the visible sector can also dilute the densities of the remaining 1 +− glueballs. We investigate these effects here, dividing the analysis into T th < T c and T th > T c cases.
T th < T c : Recall that this case is only realized for dimension-6 transfer operators, and implies that the 0 ++ decay rate is larger than Hubble following confinement. This means the 0 ++ density tracks its equilibrium value with effective temperature T x = T , and there is no longer a separately conserved entropy in the dark sector. The evolution of the 1 +− number density in this context is where n 1 denotes the equilibrium density of the 1 +− at temperature T . Note that Eq. (57) assumes the 1 +− mode also thermalizes with the visible sector. This is expected prior to freezeout since the equilibrium density of the 1 +− is smaller than that of the 0 ++ , and elastic scattering between these two species is at least as efficient as the annihilation reaction.
T th > T c : This second case implies T x ≤ T at confinement, with 0 ++ decays inactive (Γ 0 < H) until T < m 0 . To compute the resulting 1 +− relic density, we treat the 0 ++ decay as instantaneous and match the density evolution immediately before and after it occurs. Prior to the decay, the dark and visible entropies are conserved independently with ratio R, and the glueball densities evolve according to Eqs. (32,33). Decays of the 0 ++ are implemented at Γ 0 = H, where the Hubble rate includes contributions from both the visible and dark energy densities. If T i < m 0 is the visible temperature prior to the decay, the visible temperature afterwards is obtained from local energy conservation, where we have neglected the exponentially subleading contribution of the 1 +− mode to the energy density. Note that T f > T i is always smaller than m 0 as well. The evolution of the 1 +− number density after the 0 ++ decays is given by Eq. (57). Since the 1 +− number density is not changed by the decays, n 1 (T f ) = n 1 (T i ) is used as the initial condition at T = T f . The interplay of glueball annihilation, transfer, and decays leads to many different qualitative behaviours. These were investigated in Ref. [47,48] for a simplified model consisting of an unstable massive bath particle and a heavier DM state. Dark glueballs provide an explicit realization of this scenario, with the 0 ++ making up the massive bath and the 1 +− acting as (metastable) dark matter. Compared to the simple model studied in Ref. [47,48], the 0 ++ massive bath particle always freezes out (or decays) before the would-be 1 +− dark matter, corresponding to the chemical or decay scenarios discussed there. A potential further behavior that we have not captured in our approximations is the late production of 1 +− glueballs through transfer reactions while T > m 1 but after 1 +− freezeout has occurred in the dark sector. We estimate that this is potentially relevant in a very limited corner of the parameter space, and will only increase the limits we find.

C. Comments on Theoretical Uncertainties
Before applying our results for dark glueball lifetimes and densities to derive cosmological and astrophysical constraints on them, it is worth taking stock of the theoretical uncertainties in our calculations. It is also useful to identify how some of these uncertainties might be reduced with improved lattice calculations.
The glueball lifetimes computed in Sec. III rely on glueball masses and transition matrix elements. Masses for G x = SU (3) have been obtained to a precision greater than 5% in Refs. [58,59], while the matrix element relevant for 0 ++ decay was determined to about 20% in Refs. [59,69]. Thus, we expect our determination of the 0 ++ decay width to be reasonably accurate. The situation is less clear for the 1 +− width, which relies on a 1 +− → 0 ++ transition matrix element that we were only able to estimate using NDA. In the absence of lattice calculations for this matrix element, we estimate that our 1 +− width is only reliable to within a factor of a few.
Turning next to the cosmological evolution of the dark gluons and glueballs, we implicitly treated their interactions as being perturbative. This is a good approximation at temperatures well above the confinement scale, but significant deviations can arise as the temperature falls to near confinement [80]. For the range of entropy ratios R due to energy transfer computed above, this implies that values of R max with T th T c are reliable, but the specific values of R min and R max for T RH ∼ T c could receive large corrections. Similarly, the glueball interactions used to compute the (3 → 2) and (2 → 2) cross sections are quite strong for N = 3. It is difficult to quantify how this affects the pre-decay glueball relic densities, but we do note that the densities typically depend roughly linearly on R and the annihilation cross sections. Our naive estimate is that the pre-decay glueball densities we find are accurate to within about an order of magnitude.

V. COSMOLOGICAL CONSTRAINTS
In the analysis above we showed that dark glueballs can have a wide range of decay rates and a variety of formation histories in the early universe. Very long-lived dark glueballs can potentially make up the cosmological dark matter. On the other hand, shorter-lived glueballs are strongly constrained by the modifications they can induce in the standard predictions for big bang nucleosynthesis (BBN) [53,54,97], the cosmic microwave background (CMB) [56,98] and the spectrum of cosmic rays [9]. We investigate the bounds from cosmology and astrophysics on dark glueballs in this section for the four decay scenarios discussed in Sec. III. Throughout the analysis, we focus on G x = SU (N = 3), and we assume reheating such that T RH M and T x RH ≥ T c . Details of how we implement the bounds from BBN, the CMB, and cosmic rays are collected in Appendix B. To the left of the diagonal dotted lines in the lower three panels, the given fixed value of R is less than R min and is inconsistent with minimal energy transfer by the connector operators for T x RH > T c .
We see from Fig. 5 that dark glueballs are strongly constrained by cosmological and astrophysical observations. When the 0 ++ is long-lived, corresponding to small m 0 /M , its relic density tends to be too large unless the entropy ratio R is much less than unity. With sufficiently small R the 0 ++ can make up all the dark matter corresponding to the white line in the left panel of Fig. 3. Such a DM candidate would be very difficult to probe, with the most promising avenues being high energy gamma rays and modifications to cosmic structure from glueball self interactions. Using large-N and NDA, the 2 → 2 self-interaction cross section of 0 ++ glueballs is [21,25] σ 2→2 /m 0 (10 cm 2 /g) 3 N 4 100 MeV This is at (or slightly above) the current limit for N ≥ 3 and m 0 ≥ 100 MeV and could have observable effects close to these values [99], but falls off very quickly with higher mass or if the 0 ++ glueball is only a small fraction of the full DM density. For larger m 0 /M ratios, the 0 ++ and 1 +− glueballs both decay quickly enough to alter BBN or the CMB or create high energy gamma rays. Not surprisingly, the bounds from glueball decays in this scenario come primarily from the 0 ++ which has a much larger relic yield prior to decay.

B. Decay Scenario 2: Dimension-8 Decays with Exact C x
Our second decay scenario has dominant dimension-8 operators with χ i = 1 and a conserved C x charge that implies χ Y = 0 and a stable 1 +− glueball. The cosmological bounds on this scenario are shown in Fig. 6 for various values of the entropy ratio R. The upper two panels have R = R min , R max respectively, and the lower three panels show R = 10 −9 , 10 −6 , 10 −3 . As above, the grey shaded regions indicate where our theoretical assumptions are not satisfied, and the diagonal dashed lines have R < R min to their left.
The cosmological exclusions on this scenario are nearly identical to those on scenario 1 except for the new bounds from the 1 +− relic density. At the lower edge of the cyan excluded region, the 1 +− glueball can make up all the dark matter. This occurs primarily when the 0 ++ decays relatively quickly, since otherwise it tends to dilute the 1 +− relic density too strongly. Note as well The third decay scenario 3 has both dimension-6 and dimension-8 operators with y ef f = 1 and χ i = χ Y = 1, and broken C x . This leads to 0 ++ decays dominated by the dimension-6 operator, but decays of the 1 +− only through the dimension-8 operators. As a result, the 1 +− glueball is parametrically long-lived relative to the 0 ++ (and the other glueball states).
We show the cosmological and astrophysical bounds on this scenario in Fig. 7 for various values of the entropy ratio R. The upper two panels have R = R min , R max respectively, and the lower three panels show R = 10 −9 , 10 −6 , 10 −3 . As above, the grey shaded regions indicate where our theoretical assumptions are not satisfied, and the diagonal dashed lines have R < R min to their left, except in the R = R max panel. Here, thermalization is maintained all the way to confinement (and beyond) to the left of the line.
Decays of both the 0 ++ and 1 +− glueballs lead to relevant exclusions in this scenario. The 0 ++ relic density tends to be much larger than the 1 +− prior to decay, and produces the strongest constraints for small values of m 0 /M when it is long-lived. For very long lifetimes and small R, it can make up all the DM as before. However, larger values of m 0 /M lead to relatively short-lived 0 ++ glueballs that decay before the start of BBN. In this case, the longer-lived 1 +− can decay late enough to disrupt nucleosynthesis or the CMB in an unacceptable way. Note as well that the region in which the 1 +− relic density is potentially large, it decays too quickly to make up the dark matter. Our final decay scenario 4 has has both dimension-6 and dimension-8 operators with y ef f = 1 and χ i = 1, together with conserved C x (and χ Y = 0). The 0 ++ mode decays as in the previous scenario, but now the 1 +− is stable.
The cosmological bounds on this scenario are shown in Fig. 8 for various values of the entropy ratio R. The upper two panels have R = R min , R max respectively, and the lower three panels show R = 10 −9 , 10 −6 , 10 −3 . As above, the grey shaded regions indicate where our theoretical assumptions are not satisfied, and the diagonal dashed lines have R < R min to their left, except in the R = R max panel. Here, thermalization is maintained all the way to confinement (and beyond) to the left of the line.
The exclusions on this scenario from the 0 ++ are identical to those on scenario 3. However, the constraints from the 1 +− are now from its relic density rather than the effects of its decays on BBN and the CMB. This state can make up the dark matter for a range of values of its mass and the entropy ratio R. Compared to the analogous scenario 2, the relic density of the 1 +− tends to be larger here because it experiences less dilution from the more rapid decay of the 0 ++ .

VI. CONCLUSIONS
In this work we have investigated the cosmological constraints on non-Abelian dark forces with connector operators to the SM. We have focused on the minimal realization of such a dark force in the form of a pure Yang-Mills theory. In the early universe, the dark gluons of such theories confine to form a set of dark glueballs. Connector operators allow the transfer of energy between the visible (SM) and dark sectors, modify the freezeout dynamics of the glueballs, and induce some or all of the dark glueballs to decay. Late decays of glueballs can modify the standard predictions for BBN, the CMB, and cosmic ray spectra, while very long-lived or stable glueballs must not produce too much dark matter. Using these considerations, we have derived strong constraints on the existence of new non-Abelian dark forces.
A significant new feature of our work compared to previous studies [14,19,21,22,24,25,[29][30][31]44] is the inclusion of the heavier 1 +− glueball species. This state can be parametrically long-lived or stable relative to the other glueballs. It freezes out in conjunction with the 0 ++ , with the 0 ++ density forming a massive thermal bath, leading to a rich array of freezeout and decay dynamics [47,48]. In general, the (pre-decay) relic density of the 1 +− mode is much smaller than the 0 ++ . Even so, the 1 +− can sometimes yield the strongest cosmological bounds due to its longer lifetime. Specifically, the 0 ++ could decay before impacting standard cosmological processes such as BBN, while the 1 +− decays late enough to directly interfere. In some cases, the 1 +− glueball could even make up the observed DM density.
Our study also concentrated on the dark gauge group G x = SU (N = 3) with a lightest 0 ++ glueball mass above m 0 ≥ 100 MeV. The constraints found here could also be generalized to other dark gauge groups and lower masses. A very similar glueball spectrum is expected for SU (N > 3) [61], but the confining phase transition will be more strongly first-order and its effect on glueball freezeout deserves further study [24,79]. For G x = SU (2), SO(2N +1), Sp(2N ) there are no C x -odd glueballs [39,41], but otherwise we expect the constraints based on the 0 ++ glueballs to be applicable here. In the case of SO(2N > 6), the C x -odd states are expected to be significantly heavier than the 0 ++ , and thus the additional constraints on the lightest C x -odd mode would typically be weakened. Our focus on m 0 > 100 MeV was motivated by the ranges of masses considered in studies of the effects of late decays on BBN, the CMB, and gamma rays. Limits from the CMB can still be applied to much lower masses, but those from BBN and gamma ray production will be very different. We leave a study of lower glueball masses to a future work.

Cross Sections for Dimension-8 Operators
The relevant operator has the general form where F C µν is a SM field strength. This operator generates XX → AA transfer reactions for T m A , m 0 , as well as XA → XA elastic scattering. Concentrating on XX → AA, the corresponding matrix element for (p 1 , a) + (p 2 , b) → (p 3 , C) + (p 4 , D) is where a, b, C, D are "colours" and i are polarization vectors. From this expression, we find (neglecting possible masses) whereg x andg A = 2(N 2 A − 1) are the dark and visible numbers of degrees of freedom. The energy-transfer collision term is then The integral is dominated by x = √ s/T ∼ 10, corresponding to scattering at the high end of the thermal distribution.
The coupling A can be obtained by matching to our previous results for dark gluon connector operators. While there are several operators that can contribute, we keep only the S component corresponding to the operator listed above, which yields with A i = Y, 2, 3 for each of the SM gauge factors. For χ i → 1, the gluon contribution dominates withg A = 2(N 2 c − 1), and we focus on it exclusively. Note that since we are considering T m 0 100 MeV and the integration is dominated by √ s ∼ 10T , we should always be safely above the QCD confinement scale.

Cross Sections for Dimension-6 Operators
The operator of interest is now with To treat scattering through this operator, we should distinguish between temperatures above and below the electroweak phase transition at T EW P T 100 GeV. Above the transition, all the SM states are massless and we can treat the Higgs field as a complex scalar doublet. Below the transition, we must account for masses. For T > T EW P T , the dominant transfer reaction is X + X → H + H † , for which the scattering kernel is This yields the collision term ∆C =g where now the integral is dominated by √ s ∼ 6T . Below the transition temperature, we have ff , hh, W + W − , and ZZ final states at leading order. Their contributions to the scattering kernels are These results only apply for √ s > 2m i ; they are zero otherwise. Note that for √ s 2m h , 2m f , the sum of these kernels is equal to the result of Eq. (A11). electromagnetic decay products emitted before this thermalize with the photon-electron plasma before can they can destroy light elements by photodissociation [102][103][104].
The combined effects of hadronic and electromagnetic decays on BBN have been studied in a number of works, including Refs. [53,54,97]. We apply the exclusions derived in Ref. [97] to place limits on decaying glueballs, using a simple interpolation to generalize their results to arbitrary relic mass values between the range 30 GeV ≤ m x ≤ 10 6 GeV they studied, and matching to the appropriate set of final states. For masses outside these ranges, we apply the constraint for the nearest mass boundary.

Decay Constraints from the CMB
Particle decays during or after recombination at t 1.2 × 10 13 s can modify the temperature and polarization spectra of the CMB. They do so by injecting energy that increases the ionization fraction and temperature of the cosmological plasma. In turn, this broadens the last scattering surface and alters the correlations among the temperature and polarization fluctuations [56].
Detailed studies of the impact of such energy injection on the CMB have been performed in Refs. [98,[105][106][107][108][109][110]. Corresponding limits on particle decays based on the CMB measurements of Planck [84] were extracted in Refs. [98,110]. Given the theoretical uncertainties in our calculation of the pre-decay glueball yields, we apply a very simple parametrization of the results of Ref. [98]: where F(τ ) accounts for the effects of early decays. It is obtained by fitting to the curve of Fig. 4 of Ref. [98], and is normalized to unity for τ 1.2 × 10 13 s. The form of Eq. (B1) neglects mild dependences on the mass of the decaying glueball and the specific final state, but these effects are smaller than the uncertainties in the calculation of the pre-decay yield. We also apply this limit to relic masses well above the largest value of m x ∼ 10 TeV studied in Ref. [98] (and elsewhere). Such large masses lead to injections of highly energetic photons and electrons that deposit their energy very efficiently in the cosmological plasma [107]. As a result, we do not expect any major loss of sensitivity for glueball masses well above 10 TeV.
Bounds on glueball decays can also be obtained from their effects on the CMB frequency spectrum [111,112]. We find that these are subleading compared to those derived from BBN and the CMB power spectra.

Decay Constraints from Gamma Rays
Glueballs with lifetimes greater than the age of the universe t 0 4.3 × 10 17 s can produce observable signals in gamma ray and cosmic ray telescopes, even if their density is only a small fraction of the total DM value. Limits on the lifetimes of decaying DM were derived in Ref. [9] for dimension-6 glueball decays and other final states over a broad range of masses using galactic gamma ray data from Fermi [113]. With the theoretical uncertainty on glueball yields in mind, we use a simple parametrization of the limits on the glueball yield: m i Y i < (4.32 × 10 −10 GeV) τ 5 × 10 27 s e t 0 /τ e (10 GeV/m i ) , where the last two factors account for the depletion of the signal if the decay occurs before the present time and the loss of sensitivity of Fermi at lower masses [114]. This limit is fairly conser-