Indirect Detection of Dark Matter Annihilating into Dark Glueballs

We examine indirect detection of dark matter that annihilates into dark glueballs, which in turn decay into the Standard Model via a range of portals. This arises if the dark matter candidate couples to a confining gauge force without light flavours, representative of many possible complex dark sectors. Such Hidden Valley scenarios are being increasingly considered due to non-detection of minimal models as well as theoretical motivations such as the Twin Higgs solution to the little hierarchy problem. Study of dark glueballs in indirect detection has previously been hampered by the difficulty of modeling their production in dark showers. We use the recent GlueShower code to produce the first constraints on dark matter annihilating via dark glueballs into the Standard Model across photon, antiproton, and positron channels. We also fit the Galactic Centre Excess and use this observation, combined with other astrophysical constraints, to show how multi-channel observations can constrain UV and IR details of the theory, namely the exact decay portal and hadronization behaviour respectively. This provides unique complementary discovery and diagnostic potential to Hidden Valley searches at colliders. It is interesting to note that thermal WIMPs annihilating to $\mathcal{O}(10~\mathrm{GeV})$ dark glueballs and then the Standard Model via the Twin-Higgs-like decay portal can account for the Galactic Centre Excess while respecting other constraints.

1 Introduction The nature of dark matter (DM) is one of the biggest ongoing mysteries in the current era of particle physics.While decades of experiments have begun to whittle away at the parameter space of traditional WIMP-like DM, a lack of any signal of new physics across collider, direct detection and indirect detection searches has motivated interest in studying a wider variety of possible models.Inspired by the complexity we observe in our own Standard Model (SM) and the large abundance of DM relative to the SM, significant effort has recently focused on the possibility of complex dark sectors in which dark matter and other states may interact via new dark forces.
In general, the term 'dark sectors' refers to the case when we consider additional fields that can couple significantly to each other, but are uncharged under the SM.One such framework of dark sector theories are Hidden Valley (HV) models [1].Generically, the most common example includes SM-singlet fields charged under a new dark SU (N c ), see e.g.[2][3][4][5].Specific models such as Mirror Twin Higgs [6], Fraternal Twin Higgs [7], Folded SUSY [8], and many more all fall into this bracket [9][10][11][12][13][14]. Interactions with the SM are possible through portal couplings [15][16][17] that are very weak or highly suppressed at low energies.A class of dark sectors called Neutral Naturalness models, which are related to the SM via a discrete symmetry and include the Mirror Twin Higgs, can solve the Little Hierarchy Problem while evading conventional LHC searches.This makes the study of dark sectors with qualitative similarities to the complexity of the SM particularly motivated.
When dark quarks or gluons in these confining dark sectors are produced in high energy events, they will have to shower and eventually hadronize, analogous to SM QCD.This will lead to high multiplicity events equivalent to jets we observe at colliders.We refer to these events as 'dark showers' [18][19][20].A range of possible phenomena that arise from dark showers have been studied in the context of colliders, such as semi-visible [21][22][23] and emerging jets [24][25][26].Dark showers have also been studied in the context of indirect detection, such as how they explain the Galactic Centre Excess (GCE) [27,28].The effect of how phenomenologically similar cascade decays change constraints for DM searches has also been studied [29,30].While not a confining model, the phenomena of high multiplicity softened events is related.
In the N f = 0 limit, these dark showers will only contain a variety of dark glueball states [31][32][33][34][35].This is an important region of HV parameter space as it commonly appears in neutral naturalness models, but has so far been almost entirely unstudied due to the inapplicability of SM-like hadronization models [36][37][38] in the N f = 0 limit.In this work we extend the study of DM annihilating to dark showers to include this possibility, making use of recent advances in simulating pure-glue hadronization via the GlueShower code [39].
Previously dark glueballs have been considered as dark matter candidates themselves [40][41][42][43][44][45][46][47][48], but here we only focus on the case that the DM couples to the dark glueball sector, and how this affects the observed annihilation spectra compared to the traditional WIMP expectation.Additionally, we are interested to learn if these indirect detection observables allow us to glean further detailed information about the dark sector.Dark glueballs are generically long lived particles (LLPs), decaying to the SM on macro length scales.This means if they are produced in colliders they tend to decay outside the detectors.While it is possible to detect the shorter-lived states in LLP searches at the LHC main detectors [23,26,[49][50][51] and proposed external LLP detectors [52][53][54][55][56][57], states with extremely long lifetimes may be indetectable due to low production rates, or only show up in missing energy or monojet searches [58][59][60], meaning information on the decays themselves is lost.Indirect detection searches probe astrophysical length scales and could thus see the entire spectrum of glueball decays, revealing important information about the UV completion of the dark sector and/or the IR nature of pure glue hadronization.
In this study we compute the observed photon, antiproton and positron spectra in cosmic ray observations for three possible decay portals of glueballs to the SM (Higgs portal, gauge portal, and a Twin-Higgs-like mixture) and several possible benchmarks for the unknown dark hadronization physics.This allows us to provide the first constraints on the annihilation cross section in the parameter space of dark matter and dark glueball mass, and show how this simplified model can fit the Galactic Centre Excess (GCE) [61][62][63][64][65][66].We also demonstrate how multi-channel cosmic ray observations can constrain both the unknown details of pure glue hadronization and details of the dark sector UV completion by determining the decay portal for the dark glueballs.In particular, under the assumption that dark matter annihilation explains the GCE, decays via the gauge portal are excluded, while the Higgs portal and Twin-Higgs-like case are consistent, especially for O(10 GeV) glueball masses favoured in Fraternal Twin Higgs scenarios [7,67].
Our paper is outlined as follows.In Section. 2 we outline our simplified DM model we utilise in this work.We then go on to explain the various components of the dark glueball shower in Section.3, such as the properties and decays of dark glueballs, as well as how we handle pure glue hadronization.Section. 4 reviews the various methods we implement to compare our theory to data and put constraints on the DM parameter space.Lastly we provide our results in Sec. 5, the first indirection constraints for dark glueball showers, as well as Sec.6, where we outline how our analysis begins to constrains the physics of the dark glueball sector itself.We conclude in Sec. 7.

Simplified Dark Sector Model
We consider the following simplified model Lagrangian to represent the generic case of dark matter annihilating to dark gluons which then form dark glueballs: We have introduced a fermion, χ, our DM candidate; a dark confining sector containing the dark gluon G D and a dark quark q D , with mass larger than the confinement scale, Λ; and lastly a scalar ϕ that couples χ to the dark gluons via a q D loop.Minimally a single dark quark is needed for glueballs to decay to the SM, but multiple fields with a range of charge assignments and masses is possible, which we address in Section 3.2.A scalar mediator is introduced to couple the DM to the confining sector as for di-gluon production a vector mediator is forbidden by the Landau-Yang theorem for colour-singlet states [68,69].Note that a vector mediator would still allow dark glueball showers via a three gluon initial state or quirk production.We focus on di-gluon production via a scalar mediator for simplicity.Thus χ can annihilate in the galaxy to two dark gluons that then shower and hadronise into dark glueballs.Decays of dark glueballs to SM states can be allowed by coupling the dark quarks to the SM Higgs ("Higgs portal") or making them charged under SM gauge groups ("gauge portal").We discuss this in more detail in Section 3.1.While we make no particular assumption about the DM production mechanism, a minimal scenario is that χ is a thermal relic that annihilates to dark gluons in the early universe and makes up all of dark matter today.As we discuss in Section 5, it is straightforward to reinterpret our results in the more general scenario where χ is a thermal relic that has other annihilation channels and constitutes only a fraction of DM.This is important since realistic dark sectors can easily feature several types of dark relics, e.g.thermal relic abundances of stable glueball states [47] (though this never overcloses the universe for any parameters we consider and is negligible if the dark sector is cold relative to the standard model [70,71]), or asymmetric abundances of stable dark quarks in the early universe [72,73].
Our model is very theory agnostic, only containing the fields needed to produce dark glueball showers.However, it is very much inspired by the Fraternal Twin Higgs (FTH) framework [7], in which case our DM candidate χ would be the mirror tau [43,74], G D µν the mirror gluons, q D the mirror top, and ϕ the mirror twin Higgs, or, qualitatively, the twin Z/γ. 1 While the original twin tau DM proposal has since been excluded by direct detection constraints, it was recently understood [75] that including the effects of spontaneous Z 2 breaking in the FTH setup makes thermal twin tau dark matter a viable candidate and an attractive target for future direct detection experiments.
Note that in the case of q D coupling to both ϕ and h, scalar mixing will allow χ to annihilate directly to the SM.However we are interested in the case when this coupling or mixing is small, thus this channel will be suppressed compared to dark gluon production.Small y q values would also increase the lifetime of the dark glueballs.

Production and Decays of Dark Sector Glueballs
In this section we briefly summarise the known physics of dark glueballs.Firstly, we discuss the properties of the various glueball states, such as masses and quantum numbers.We then explain how the dark glueballs are able to decay to the SM via higher dimension operators and how we calculate their branching ratios.Lastly we outline what we know about the production of dark glueballs via pure glue hadronization, its implementation in GlueShower, and importantly how we address uncertainties due to the unknown pure-glue non-perturbative physics.

Basics
The properties of SU (N c ) glueballs have been studied on the lattice for decades [31][32][33][34][35], establishing a spectrum of twelve stable states in the absence of external couplings, as shown in Figure 1.The spectrum of masses is given in terms of m 0 , the lightest glueball mass, and labelled by their quantum numbers J P C .Alternatively the masses can be parameterised by the confinement scale of the theory, Λ, since it has been calculated on the lattice that m 0 ≈ 6Λ.Thus there is only one free parameter, m 0 , in determining the entire mass spectrum of glueballs.In this work we will only study the case N c = 3 2 and use the 1 To be precise, in the FTH model the annihilation of twin tau DM to mirror gluons occurs via twin higgs and Z bosons in the s-channel and a dark quark loop to either two or three gluons directly, or to an intermediate quirky [2] bound state of mirror bottoms, which then annihilates to mirror gluons.The details of quirky de-excitation and annihilation in this case may be complicated, but the outcome for glueball production should be very similar to the simplified scenario we study, since glueball production should be dominated by s-wave annihilation of the mirror bottoms.
2 Considering different values of Nc would only slightly change the mass ratios for the various glueball states and not significantly change our results.
Figure 1.Glueball mass m G spectrum for pure SU (3) Yang-Mills theory [31] in terms of the lightest glueball mass m 0 , plot taken from [76].
mass values calculated in Ref. [34], only considering the lightest 10 glueballs.We focus on m 0 values in the range 10 -50 GeV, as this regime is motivated by theories of neutral naturalness [7,67], but will consider an extended range of glueball masses for completeness.Depending on the charge assignment of the dark quark, there are two possible sets of effective operators that couple the dark sector to the SM [76,77].Firstly, if the dark quark couples to the Higgs field, we get an effective dimension-6 operator: where M is the mass scale of the dark quarks, H the SM Higgs field, and G D the dark gluon field strength.In the Twin Higgs model, this arises due to a small required mass mixing between the twin and visible Higgs mass eigenstate.Secondly, if the dark quark is charged under any of the SM gauge groups, we also get dimension-8 operators: where G SM is a SM field strength.In the Twin Higgs model, these operators may arise from a UV completion which often features fields above a TeV charged under both sectors [78][79][80].These operators allow some or all of the dark glueballs to decay to the SM.Since these operators can be factorised into respective SM and dark sector operators, we need to consider the forms: where D is the overall dimension of the operator and Θ κ refers to a dark glueball with quantum numbers labelled by κ.The SM part ⟨SM |O (D−d) |0⟩ can be evaluated using perturbative field theory methods, but the dark sector transition matrix elements require non-perturbative physics.We will discuss the specifics of the required matrix elements and decay channels that arise from these operators in the following subsection.

Decay Modes
The relative strengths of the aforementioned Higgs portal and gauge portal operators are effectively a free parameter of our model.We analyze three particularly simple or motivated cases: • Higgs Portal: In this regime the dark quark only couples to the Higgs, thus decays only proceed through the dimension-6 effective operator.The majority of the dark glueballs are able to decay to the standard model, this is achieved by either directly mixing with the SM Higgs or by decaying to a lighter glueball and radiating an off shell SM Higgs.Note that the 0 −+ and 1 +− states are entirely stable.
• Gauge Portal: In this regime the dark quark only couples to the SM gauge fields, thus decays only proceed through the dimension-8 effective operator.The lighter glueballs decay to two SM gauge bosons, while the heavier states mostly decay to a lighter glueball and radiate a single photon, or a Z boson if allowed.In this regime none of the dark glueballs are stable.
• Twin Higgs-like: In this regime we assume both operators are active but we introduce two different scales, M (6) , the mass scale of the dimension-6 operator, and M (8) , the mass scale of the dimension-8 operator.We consider the limit M (6) ≪ M (8) , inspired by UV completions of Twin Higgs theories [78][79][80].A mirror top partner, whose mass sets M (6) , is present and couples to the SM Higgs but not the SM gauge bosons.At some higher scale, M (8) , the theory is UV completed, with fields that have SM gauge charge.In this case due to the mass hierarchy and relative dimensions of the effective operators, dimension-6 decays will heavily dominate over dimension-8 decays.Thus when both decay channels are allowed for a glueball species, we only consider dimension-6 decays, but the previously stable 0 −+ and 1 +− states will now decay via the dimension-8 operator.
We summarise the decay modes in Table 1.The specific formula for the various decays can be found in Refs.[76,77] but we will now discuss some of the most important features of this decay phenomenology.
For the Higgs portal case, the 0 ++ glueball mixes directly with the Higgs and has branching ratios: where ++ ⟩ is the 0 ++ annihilation matrix element, and Γ SM h→ξξ (m 0 ) is the SM Higgs branching ratios as a function of mass, calculated using HDECAY [81] or chiral perturbation theory [82] (see below).Thus for the 0 ++ state, the branching ratios are equivalent to a SM-like higgs of mass m 0 .However, at larger m 0 values, the 0 ++ → hh channel also becomes available.
The rest of the Higgs portal decay modes take the general form, Θ κ → Θ κ ′ +h * , with the off shell Higgs then decaying to the SM.To calculate the branching ratios for these decays, we need the matrix transition elements, ⟨Θ κ ′ |tr(G D G D )|Θ κ ⟩.While some of these matrix elements have been calculated on the lattice [83,84], the majority of them remain unknown.However the various matrix elements these studies calculate are within O(1) factor of one another, so we make the approximation that for the purposes of branching ratios, these matrix elements factorise out and cancel.Under this assumption, the branching ratios are determined solely by the mass splittings between the glueball states.In the m 0 range where h * is produced far off shell, we treat the decay in two steps.Firstly, by binning the two body phase space as a function of m h * , we are able to define branching ratios to h * of various mass values.Once the glueball has decayed to the off shell Higgs with a specific mass value, we use HDECAY (if m h * > 2 GeV) or chiral perturbation theory [82] (if m h * < 2 GeV) to determine its branching ratios into SM fundamental particles or hadrons.In the case that the Higgs can be produced on shell, it is treated as a simple two body decay.
For dimension-8 operators, the branching ratios to different gauge groups is dependent on the representation(s) of the heavy field(s).This is encoded in the various coefficients that appear in the full Lagrangian [47,76,77]: where r is the SM representation of the heavy dark quarks with masses M r , and ρ r = M r /M .Index, i, runs over the sub-representations corresponding to the SM gauge groups: indicates the number of copies of the i-th subrepresentation, and T 2 (r i ) is the trace invariant for that factor.To remain relatively theory agnostic about the UV completion of the model, we simply set χ i = χ Y = 1.This gives us the relation: which applies for Θ = 0 ++ , 2 ++ , 0 −+ , 2 −+ .For small m 0 values, these are the only significant decay channels for the dark glueballs listed.Thus the dimension-8 branching ratios for the lighter glueballs are entirely determined by SM parameters.As m 0 increases, the Zγ, ZZ, W W channels become available and their branching ratios can be found in Ref. [76,77].
The majority of the heavier glueball dimension-8 decays are of the form, Θ κ → Θ κ ′ + {γ, Z}.Once we assume the corresponding matrix elements are of the same order and cancel, the branching ratios are determined by the glueball mass splittings and a spin dependent factor.
The 3 ++ state is noteworthy as it has three competing decay modes due to the dimension 8 operator, and a large uncertainty regarding which dominate.In this work we simply assume 3 ++ → 1 +− + {γ, Z}, ignoring the 3 ++ → {2 ++ , 0 −+ } + {γ, Z} decay modes.We select this decay mode as it is the most distinct and will provide the most conservative constraints.Additionally this will be a small effect on the total spectrum.The 3 ++ state is a heavy glueball and produced in relatively small amounts compared to the lighter glueballs across the hadronization benchmarks, as we explain in the following subsection.

Pure Glue Hadronization
To calculate the energy spectra of the SM particles produced in dark glueball showers, we first need the energy spectra of the various dark glueball states produced in DM annihilation.These spectra, also known as fragmention functions, are dependent on the model used to shower and hadronize the initial dark gluons into the final state dark glueballs.In this work we will use GlueShower3 [39] to generate these dark glueball fragmentation functions.
We will briefly summarise the underlying physics GlueShower utilises in its algorithm, see Ref. [39] for a more in-depth discussion.Firstly, like any SU (N c ) theory, the shower evolution of the initial virtual dark gluons is determined by perturbative QCD.This allows us to calculate the virtualities and multiplicity of dark gluons as the shower evolves towards the confinement scale.Assuming that N f = 0 QCD satisfies the assumption that hadronization is an IR process active around the confinement scale (referred to as the 'jetlike' case), the shower is then terminated not at the confinement scale but at Λ had ≥ 2m 0 .This is due to the fact perturbative estimates of momentum exchange between separate branches of the shower past this point were found to be very small, and therefore such a dark gluon can no longer split and form two on shell glueballs once confinement occurs.The dark gluon is hence simply put on shell as a single dark glueball.

Hadronization Hadronization Hadronization Shower
In both cases described above, the relative probabilities of the various glueball states are determined by a Boltzmann distribution, controlled by the hadronization temperature, T had , calculated from lattice QCD [85].Thus the production of heavier glueball states are exponentially suppressed.
In addition to the qualitative choice of shower behaviour, two nuisance parameters are included to account for the unknown quantitative details of pure-glue hadronization: Λ had /2m 0 ≥ 1 affects the final state multiplicity of dark glueballs, and T had /T c ≥ 1 determines the relative production of heavier vs lighter dark glueball states.The dependence of glueball fragmentation functions on these nuisance parameters is explored in Ref. [39] and they provide a set of eight benchmark parameter choices that bracket the possible range of phenomena.In Table .2 we list the eight benchmark sets of nuissance parameters and an example of the resultant fragmentation functions corresponding to these benchmarks are shown in Fig. 2. The range in behaviour of these fragmentation functions is used in our study to incorporate the theoretical hadronization uncertainties into our final constraints.
An interesting limit to consider is when the centre of mass energy of the shower E is less than the scale at which we end the perturbative shower, 2m 0 ≤ E ≤ Λ had .For the jetlike case, each gluon is not allowed to perturbatively split as its virtuality is already below the hadronization scale.Thus for E ≤ Λ had , jet-like showers always result in back-to-back glueball production.For the plasma-like case, again the gluons are unable to perturbatively split.However, in this case the gluons are then treated as a single plasma at rest that deexcites quasi-isotropically as mentioned before.

Methodology
To study DM annihilating into dark glueballs, we first must calculate the resulting SM spectra in the galactic rest frame.We show this for a variety of dark matter and glueball masses, as well as across our three decay portals and eight hadronization benchmarks.Once an assumption is made regarding the distribution of DM in our galaxy, we can then calculate the observed flux at Earth's position.For charged particles this requires cosmic ray propagation and solar modulation tools.Finally with the flux at Earth calculated, we are able to find constraints by comparing the spectra to observational data by Fermi-LAT in the case of photons, and AMS-02 for antiprotons and positrons.

Monte Carlo Generation of SM final state spectra in DM rest frame
For each species of dark glueball, we generate the energy spectra of photons, positrons, and antiprotons from their decays using PYTHIA 8 [86].The fragmentation function for each dark glueball species, produced in the annihilation of DM to dark gluons, is generated using GlueShower.We then Lorentz boost and convolve the SM spectra with the dark glueball fragmentation functions to get the SM energy spectra for the entire dark glueball shower in the galactic frame.
Fig. 3 depicts the photon spectra in the galactic frame for a reference dark matter and glueball mass.In the simplest Higgs portal case we see how our assumptions about hadronization affect the spectral energy peak position as well as total flux.Plasma-like showers lead to softer spectra compared to jet-like showers, and how higher temperatures lead to overall smaller multiplicity due to an increased production of heavier, stable glueballs.In both the gauge portal and Twin Higgs-like cases we see the effect of heavier glueballs decaying to lighter glueballs and a single photon.Due to the range of produced glueball energies, this spectral line is smoothed out to a second, high energy peak in the spectra.The presence and height of this second peak is largely influenced by the hadronization temperature.
We also show the antiproton spectra, in the galactic frame, in Fig. 4. We find there is less of a dependence on the decay portal determining the spectral shape, compared to   the photon channel.However, there are larger differences in peak positions and total flux depending on the hadronization assumption used for the antiprotons.This is to be expected as the antiproton's energy is more influenced by the boost provided by the parent glueball due to its large mass compared to the photon.Thus the antiproton channel is more sensitive to the glueball fragmentation function and therefore hadronization assumptions.Once we have calculated the SM spectra in the galactic frame, to get the actual observed flux at Earth, we must make an assumption regarding the distribution of DM in the galaxy, which we will describe in the following section.Additionally for positrons and antiprotons, their galactic frame spectra have to then be propagated to the solar system's position using a cosmic ray propagation tool.This is due to their trajectories being influenced by the galactic magnetic fiields.Commonly used examples are fully numerical codes such as GALPROP [87], DRAGON [88,89], and PICARD [90], or even semi-analytical tools such as USINE [91].Once the positron and antiproton spectra have been calculated at the solar system's position, we then consider solar modulation effects from the heliosphere.We outline our method of cosmic ray propagation and solar modulation in Sec 4.2.

Cosmic Ray Propagation
The evolution of the number density spectra, N i (⃗ r, p), throughout the galaxy for each particle species we consider is given by the diffusion equation [92,93], where D(R) and D pp (R) are the spatial and momentum diffusion coefficients respectively.Apart from spatial diffusion and momentum diffusion (also referred to as reacceleration), this equation also includes the effects of convection, energy losses, injection from sources, and decays and collisions with the interstellar medium.
In this work we use a customised version of the DRAGON2 code4 to solve the CR diffusion equation and calculate propagated energy spectra for various particle species.This code uses a spatial diffusion coefficient with a rigidity dependent break, taking the form where R 0 = 4 GV .Recent studies have shown that CR data favours a diffusion coefficient with a high-rigidity break [94], and from this work we use their parameterisation with R b = 312 GV , ∆δ = 0.14, s = 0.04.
The momentum diffusion coefficient takes the form which is related to the spatial diffusion coefficient, as well as the Alfvén velocity, V A .For the remaining free parameters we implement the values found in Ref. [95], which also used the customised DRAGON2 code and found these values by performing a combined analysis, fitting multiple simulated CR ratios and spectra to observed cosmic ray data.
The DM source contribution is given by: where ρ(⃗ r) is the DM halo density profile, ⟨σv⟩ is the thermally averaged interaction cross section, and dN/dE is the energy spectra from DM annihilation.To account for solar modulation once the spectra have been propagated with DRAGON2, we use the force-field approximation [96], characterised by the Fisk potential, ϕ: However, solar modulation is known to be both time and charge dependent, and attempts have been made to use more sophisticated solar modulation techniques.To incorporate some of these effects we use a charge-dependent Fisk potential: where ϕ ± 1 is a charge-dependent term that accounts for increased energy loss when a particle's sign opposes the solar polarity.Thus we assume ϕ + 1 = 0, and let ϕ − 1 be non-zero.While this is slightly more complex than the basic force-field approach, we note this is still an approximation and there exist many more involved methods where the solar propagation is solved numerically [97][98][99].

Comparison to Data
For the dark glueball decay chains we consider in this work, the majority of SM states that are produced are hadronic, arising from the production of bb or gg pairs in dark glueball decays.Searches for antiprotons and gamma-rays produced by these events provide the best option for finding a potential signal but also put the strongest constraints on the DM parameter space.Additionally we also include constraints from positron searches for the sake of completeness and as an extra search channel, even though they provide less competitive constraints.In the following sections we outline how these constraints are calculated from the observed data.

Galactic Centre Excess
Various interpretations of the status of the GCE exist in the literature 5 , with different analyses changing the shape of the perceived excess [63][64][65].We fit our DM flux, given by: to the excess as calculated in Ref. [65].The value of J GCE is determined by both the region of interest and DM density profile.Following the analysis of Ref. [100], we assume a ROI of 40 • × 40 • centred on the galactic centre, and take their MED DM density profile parameterisation.They use a generalised Navarro-Frenk-White (NFW) function [101], of the form: 5 See Ref. [66] for a comprehensive review where γ = 1.3, ρ s = 0.449 GeV cm −3 , and r s = 12.67 kpc.This DM density profile is also used for the DM source term when using DRAGON2 to propagate the antiproton and positron spectra.

Fermi-LAT Dwarf Spheroidal Gamma-ray Limits
We follow the official Fermi analysis on Pass 8 LAT data [102] to find constraints on DM annihilating to photons via a dark glueball shower.This method is a combined analysis of 41 galaxies, 28 of which are kinematically confirmed, while the remaining 13 are considered likely galaxies.
Fermi publicly releases their bin-by-bin likelihood functions6 , as a function of energy flux for each dwarf galaxy.The uncertainty in the J-factor is also included in our likelihood function as a nuisance parameter we profile over, following the profile likelihood method [103].This modifies the likelihood for each galaxy in the following manner: where J i is the true value of J-factor, while J obs,i is the measured J-factor with error σ i .Values for J obs,i and σ i are taken from Ref. [102] when provided, but in the case that the data is not available we use the predicted value, J pred,i , with a nominal error of 0.6 dex, assuming a NFW profile.Constraints on ⟨σv⟩ for a given m DM are calculated using the delta-log-likelihood technique, which for a 95% C.L upper limit requires a change of 2.71/2 from the maximum value.

AMS-02 Antiprotons
We compare our generated spectra to the latest antiproton spectra released by AMS-02 [104].We assume the data observed by AMS-02 contains no DM source contributions, i.e. the spectra is solely SM background.While there has previously been claims of an antiproton excess that could be explained by a DM component [105,106], establishing the significance of these claims has been difficult since the error correlation matrix of the AMS-02 analysis is not publicly available.Recent analyses using the updated data we consider in this work [100,107], or implementing their own error covariance estimates [108,109], were able to reduce the significance of the signal or even make it go away.It is from these developments we make the previous assumption that the observed data is solely SM background.
From this assumption we choose to keep the DRAGON2 propagation parameters fixed for all DM parameter choices as they best fit the observed SM data.Due to the uncertainty in solar modulation and antiproton cross section, we let the solar modulation parameters and SM antiproton spectra normalisation, N , float along with the DM signal normalisation to find the best fit.From combined fits to the AMS-02 [110] and Voyager proton data [111], the Fisk potential parameter was determined to lie within the range, ϕ 0 = 0.60 − 0.72, which we let it float within.ϕ − 1 is allowed to float freely.There is roughly 10% uncertainty on the secondary p cross section [112], important to the SM background processes where p are produced by protons scattering off the interstella medium.Due to this uncertainty we let N float in the range 0.9-1.1.
Additionally, to minimise the influence of assumptions made about solar modulation, which is known to only have large effects on low-energy cosmic rays, we only fit the spectra above 4 GeV.The best fit is determined by minimising: where θ is the set of our theory parameters.The DM signal normalisation is then increased until the fit worsens such that the χ 2 value increases by 2.71, corresponding to a 95% C.L. upper limit on ⟨σv⟩.

AMS-02 Positrons
Similar in process to the antiproton constraints, we use the latest AMS-02 positron data [113] and assume it contains no DM contributions.However, for the positron spectra AMS-02 provides an analytical function that provides a good fit to the data.We use this function as our SM background and add a DM contribution found using DRAGON2.While we fit the DM contribution to the observed data, the analytic function's free parameters are allowed to float within their significance range.Again the quality of the fit is determined by χ 2 (θ) and a 95% C.L. upper limit on ⟨σv⟩ is found by increasing the DM signal normalisation until χ 2 (θ) increases by 2.71.

Constraints on Dark Matter annihilating to Dark Glueball Showers
Here we present the 95% confidence limits on DM ⟨σv⟩ across photon, antiproton, and positron channels.Photon constraints are found by comparing to Fermi-LAT dwarf spheroidal galaxy observations, while both antiproton and positron constraints are found by comparing to AMS-02 cosmic ray data.To incorporate the theoretical uncertainty of our hadronization models used, we calculate the confidence limit for each of the eight hadronization benchmarks.Treating the jet-like and plasma-like cases separately, due to their largely different underlying assumptions, we then depict a single constraint limit for each case, with a width defined by the maximum and minimum values across the respective four hadronization benchmarks we consider.
In the following plots we include the thermal relic cross section ⟨σv⟩ thermal [114] as a benchmark.We show excluded annihilation cross sections, ⟨σv⟩ limit , which are computed assuming that χ makes up all of the DM and annihilates exclusively into dark gluons, but our results can be interpreted more generally.If the sensitivity exceeds the thermal benchmark ⟨σv⟩ limit < ⟨σv⟩ thermal , then this effectively constraints showers with m 0 = 10 GeV (top) and 50 GeV (bottom).We plot the limits for jet-like (blue) and plasma-like (orange) showers separately for each of the decay portals we consider.The thermal relic cross section is the black dashed line [114].
in the more general scenario where χ may only constitute some mass fraction of DM and may only produce dark glueballs in a fraction of its annihilations.On the other hand, ⟨σv⟩ limit > ⟨σv⟩ thermal for a given set of parameters means we are not sensitive to any thermal relic χ regardless of its DM mass fraction or annihilation fraction to dark glueballs.

Fermi-LAT Dwarf Spheroidal Galaxies
Figure 5 shows the 95% confidence level limits found by comparing our DM photon spectra to dwarf spheroidal galaxy data.Here we focus on m 0 values 10 and 50 GeV as they bracket the motivated rage for neutral naturalness theories.For m 0 = 10 GeV, across all decay portals we consider, jet-like showers are more constrained compared to the plasma-like case.Thus even with plasma-like showers resulting in higher multiplicity glueball events, the softer spectra leads to weaker constraints.Comparing decay portal cases, the presence of an active dimension-8 operator also leads to stronger constraints.This is to be expected as this operator leads to the production of high energy gamma ray emission, compared to the softer spectra produced in hadronic showers present in dimension-6 operator decays.
Looking at the behaviour of the confidence limits themselves, as expected we also see the jet-like and plasma-like constraints converge to the same limit as m DM approaches m 0 .This is because in this limit, regardless of hadronization assumption, by energy conservation the Figure 6.AMS-02 antiproton limits at 95% C.L. for DM annihilation to dark glueball showers with m 0 = 10 GeV (top) and 50 GeV (bottom).We plot the limits for jet-like and plasmalike (orange) showers separately for each of the decay portals we consider.The thermal relic cross section is the black dashed line [114].
shower can only produce two 0 ++ glueballs.As m DM begins to increase, initially the shower limits remain similar across benchmarks, as they still can only produce two glueballs, but now producing the heavier states.However above m DM = 2m 0 the jet-like and plasmalike limits begin to diverge as hadronization assumptions start to play a role on spectra evolution.The jet-like shower limit quickly widens as hadronization dependent behaviour becomes influential.When all jet-like benchmarks have passed their hadronization scale, the jet-like limit maintains a relatively constant width and curve.At large energies the bottom of the band begins to plateau, which is a known feature of dwarf spheroidal constraints [102].For the plasma-like shower, the limit quickly weakens compared to the jet-like case.This is because the shower is still in the limit that it is hadronizing as a single plasma that then emits thermal momenta glueballs.Thus even though m DM is increasing, the glueballs are only produced with on average a constant thermal momenta, leading to weakening constraints.However, once 2m DM > Λ had , the gluons perturbatively split before plasma production, leading to boosted plasma and final state glueballs with increasing energy.In this regime you see the limit change in behaviour and approach the jet-like limit.For heavier glueballs with m 0 = 50 GeV, the constraints are slightly weakened compared to the light glueball limit for the same m DM .Due to this fact and only probing DM masses above 50 GeV, the constraints do not exclude thermal relic DM of any mass.
We also point out that even with considering the eight various hadronization benchmarks, the resultant jet-like and plasma-like constraint bands have a modest thickness and lie mostly on top of each other.Therefore, Fermi-LAT Dwarf Spheroidal limits are fairly insensitive to the unknown details of glueball hadronization.

AMS-02 Antiprotons
Figure 6 shows the 95% confidence limits found by comparing our DM antiproton spectra to AMS-02's antiproton data, again for the motivated bracketing glueball masses m 0 = 10 GeV and 50 GeV.We immediately see that these constraints are generally much stronger than the dwarf spheroidal galaxy constraints, which is expected since dark glueballs predominately produce hadronic showers over direct gamma rays.It is also clear that unlike the Dwarf Spheroidal gamma ray limits, antiproton limits are much more sensitive to the unknown aspects of dark glueball hadronization for intermediate glueball masses, and the jet-or plasma-like assumptions lead to qualitatively different behaviour.
While looking quite different than the dwarf spheroidal constraints, many aspects of the previous analysis are still applicable.Some new behaviour is that in the limit m DM approaches m 0 , the constraints significantly weaken.This is due to the fact the showers produced are so soft they are in the energy regime highly dependent on solar modulation effects.Since we only fit the spectra above 4 GeV to remain conservative about solar modulation assumptions, we lose sensitivity for this m DM range.
Again, until around m DM = 2m 0 , both limits remain relatively inline, but past this point begin to diverge drastically as the plasma-like limits severely weaken again.This feature is again due to the antiproton spectra being produced with the same average thermal momenta while the dark matter mass increases.However, once the hadronization scale is saturated, at 40 and 60 GeV respectively, the glueball energy begins to increase and the limits increase in strength, approaching the jet-like constraints in the large m DM limit.The m 0 = 50 GeV case follows the same qualitative behaviour.

AMS-02 Positrons
Figure 7 shows the 95% confidence level limits found by comparing our DM positron spectra to AMS-02's positron data, for both light and heavy dark glueballs.These limits are included for completeness, as they do not provide any competitive limits compared to both the photon and antiproton channels.A unique feature of the positron limits however is that the plasma-like showers provide stronger limits for larger DM masses.In this regime the plasma-like spectra have a larger energy flux in the low energy bins, compared to the jet-like case, thus providing stronger constraints.

Combined Constraints
In this section we summarise the constraints provided by all observation experiments and present them in the full (m DM , m 0 ) mass plane in Figure 8 for all decay portals and hadronization benchmarks.We plot the strongest constraint for each point in our parameter space, which generally is the AMS-02 antiproton limit, except for the small m 0 Higgs Portal     m 0 = 50 GeV, with best fit values residing in the range 3.45 -6.55.For each hadronization benchmark, the simulated spectra are able to fit the main peak of the GCE, but lack the high energy tail that is also present.
In comparison, the gauge and Twin Higgs-like portals are able to provide a better fit to the high energy tail, but only for the low hadronization temperature (T had = T c ) showers.For m 0 = 10 GeV, these low temperature χ 2 /d.o.f values lie in the range 4.13 -9.26 and 4.42 -9.79 respectively.The high temperature showers are unable to provide a good fit due to the presence of a large second spectral peak, causing increased χ 2 /d.o.f values in the range 68 -123 and 16.70 -37.15.For this reason we will discount them from the rest of our analysis.The m 0 = 50 GeV case is similar in the fact that the high temperature showers exhibit a second, sharp peak that does not fit the shape of the high energy tail.
In Figure 10 we plot the favoured DM parameter space in the m DM − ⟨σv⟩ plane and the previously calculated constraints for the Higgs portal case.For the jet-like showers, Twin-Higgs-like, m 0 = 10 GeV   Combined plot for the Twin Higgs-like decays, with m 0 = 10 GeV (top) and 50 GeV (bottom), containing limits at 95% C.L from both AMS-02 antiproton data (solid) and Fermi-LAT dwarf spheroidal data (hatched).The 3σ contour for the best fit to the GCE is also depicted in purple for each hadronization benchmark.The thermal relic cross section is the black dashed line [114].
only the low temperature showers provide a good fit to the GCE and are further considered.For m 0 = 10 GeV, both jet-like and plasma-like showers are still allowed, with even the GCE favoured parameter space for P 4,1 agreeing with the thermal relic cross section.For m 0 = 50 GeV, the jet-like shower begins to be constrained as J 2,1 no longer can the GCE and the antiproton limits.This trend continues for larger glueball masses with the favoured parameter space continuing to be constrained.We check glueball masses for values above the motivated range, 10 -50 GeV, and find that all low temperature showers get excluded as the spectral peak eventually shifts to larger energies and can no longer reasonably fit the GCE.
It is interesting to note that dark matter annihilating into glueballs decaying to SM particles via the Twin-Higgs-like combination of Higgs and Gauge portals can explain the GCE with a near-thermal annihilation cross section while being consistent with other constraints.Antiproton data in particular seems to favour glueball masses in the lower range near 10 GeV, which is motivated by RG arguments for Fraternal Twin Higgs models [67].This supports the compelling possibility of a WIMP-like dark matter candidate embedded in the Twin Higgs framework, for example a ∼ 50 − 100 GeV mirror tau [75].

Conclusions
Complex dark sectors are extremely motivated models that are able to provide dark matter candidates as well as mechanisms to solve the Little Hierarchy problem.A general consequence of these models is the existence of long lived particles that are being targeted by LHC searches [23,26,[49][50][51] and proposed dedicated external LLP detectors [52][53][54][55][56][57].In this work we show that cosmic ray searches for dark matter annihilating into dark glueballs can provide unique and complementary handles for both discovery and diagnosis of such dark sectors.In the regime where some or all of these long lived particles decay outside of our detectors, information about their decay is lost.Since indirect detection probes astrophysical length scales, there is a greater chance to see the entire spectrum of decays, especially if the dark sector contains a range of particles of varying lifetimes, as is the case for dark glueballs.
We utilised recent progress on simulating pure glue hadronization [39] to present the first indirect detection constraints on dark matter annihilating to dark glueballs that decay to the Standard Model, see Figure 8.We provide these constraints taking into account uncertainties arising from the unknown nature of pure-glue hadronization, and for three choices of portals via which the dark glueballs decay: Higgs portal, gauge portal and a Twin-Higgs-like mixture, corresponding to different physics in the UV.Despite the uncertainties, we derive robust constraints on the masses of dark matter and dark glueballs, with photon and antiproton constraints probing thermal annihilation cross sections.
Perhaps our most interesting finding is that multi-messenger observation of dark matter annihilation can be used to effectively constrain and even diagnose the properties of the dark sector beyond the masses of the dark matter and dark glueball.We found that the DM photon spectra are most sensitive to the dark glueball decay portals, while the DM antiproton spectra are most sensitive to hadronization behaviour.Figures 10,11 and 12 show best-fit regions of the GCE alongside Dwarf Spheroid and antiproton constraints for the Higgs portal, gauge portal and Twin-Higgs-like case respectively.Gauge portal decays seem to be entirely excluded from explaining the GCE.On the other hand, the Higgs portal and in particular the Twin-Higgs-like portal are consistent with generating the GCE for thermal annihilation cross sections and plasma-like dark hadronization behaviour, especially for O(10 GeV) glueball masses that are favoured in Fraternal Twin Higgs models.
The possibility of a thermal WIMP-like dark matter candidate embedded in a Neutral Naturalness framework is compelling for many reasons [43,74,75].Our work shows for the first time that this scenario is explicitly compatible with cosmic ray constraints and can specifically account for the GCE.If this scenario were realized by nature, then a plethora of complementary collider and astrophysical observations could allow us to detect dark matter, the dark sector as a whole and its connection to the little Hierarchy problem, an exceedingly exciting prospect.This also motivates understanding other subtle signatures of such confining dark sectors in more detail, like the enhanced production of anti-nuclei [115], which could explain the observation of eight antihelium events by AMS-02 [116], a smoking gun signal for DM [117,118].However, the more general lesson of our study is just as optimistic: complex dark sectors produce complex indirect detection signatures, with enough handles to extract many aspects of the dark sector physics that may be inaccessible at colliders, providing unique and complementary discovery and diagnostic potential.With future telescopes such as GAPS [119] launching in the near future, indirect detection searches for dark matter are about to enter a very exciting era.

= 1 ,Figure 3 .
Figure 3. Photon spectra in the galactic rest frame for m DM = 100 GeV and m 0 = 10 GeV.Depicted are the spectra for each of the eight hadronization benchmarks and three decay portal cases we consider.

Figure 4 .
Figure 4. Antiproton spectra in the galactic rest frame for m DM = 100 GeV and m 0 = 10 GeV.Depicted are the spectra for each of the eight hadronization benchmarks and three decay portal cases we consider.

Figure 5 .
Figure 5. Fermi-LAT Dwarf Spheroidal limits at 95% C.L. for DM annihilation to dark glueball showers with m 0 = 10 GeV (top) and 50 GeV (bottom).We plot the limits for jet-like (blue) and plasma-like (orange) showers separately for each of the decay portals we consider.The thermal relic cross section is the black dashed line[114].

Figure 8 .
Figure 8. Strongest indirect detection constraints for Higgs portal (top), gauge portal (middle), and Twin Higgs-like decays (bottom).For each decay portal we show constraints for each of the 8 hadronization benchmarks in Table.2.Red dotted line indicates where the strongest constraint switches from Fermi-LAT (lower left) to AMS-02 antiprotons (upper right).

Figure 12 .
Figure 12.Combined plot for the Twin Higgs-like decays, with m 0 = 10 GeV (top) and 50 GeV (bottom), containing limits at 95% C.L from both AMS-02 antiproton data (solid) and Fermi-LAT dwarf spheroidal data (hatched).The 3σ contour for the best fit to the GCE is also depicted in purple for each hadronization benchmark.The thermal relic cross section is the black dashed line[114].
Table of decay channels for each glueball we consider.h * indicates an off shell Higgs, but if allowed will be produced on shell.

Table 3 .
Best fit values for DM parameters found by fitting to the GCE spectra reported in Ref.