Gamma-rays from Dark Showers with Twin Higgs Models

We consider a twin WIMP scenario whose twin sector contains a full dark copy of the SM hadrons, where the lightest twin particles are twin pions. By analogy to the standard WIMP paradigm, the dark matter (DM) freezes out through twin electroweak interactions, and annihilates into a dark shower of light twin hadrons. These are either stable or decay predominantly to standard model (SM) photons. We show that this 'hadrosymmetric' scenario can be consistent with all applicable astrophysical, cosmological and collider constraints. In order to decay the twin hadrons before the big-bang nucleosynthesis epoch, an additional portal between the SM and twin sector is required. In most cases we find this additional mediator is within reach of either the LHC or future intensity frontier experiments. Furthermore, we conduct simulations of the dark shower and consequent photon spectra. We find that fits of these spectra to the claimed galactic center gamma-ray excess seen by Fermi-LAT non-trivially coincide with regions of parameter space that both successfully generate the observed DM abundance and exhibit minimal fine-tuning.


I. INTRODUCTION
The twin Higgs (TH) mechanism provides a natural means to stabilize the electroweak scale up to O(10 TeV) with uncolored top partners, that belong to a twin dark copy of the standard model (SM) [1,2]. This framework thus can relax the evolving tension between naturalness expectations and the current limits on traditional top partners, and is therefore an attractive solution to the little hierarchy problem.
The central idea behind the TH mechanism is that the Higgs is realized as a pseudogoldstone boson of an accidental global symmetry, which suffices to protect its mass from one loop quadratic corrections. In the simplest case, this structure is the result of an (approximate) Z 2 exchange symmetry between the SM and twin sectors. (See Refs. [3,4] for a systematic discussion of more exotic options.) In order for the Z 2 exchange symmetry to provide sufficient protection [5], the matter content of the twin sector must at the very least contain a twin third generation of quarks -a twin top and twin bottom -charged under twin color SU (3) c and twin electroweak SU (2) L gauge groups, such that the twin top yukawa and twin SU (2) L coupling are close to the corresponding SM values. Naturalness does not significantly constrain the remaining features of the twin sector, but all TH models contain a lightest twin particle (LTP) and a lightest twin hadron (LTH). The latter is either a bound state of twin quarks or a twin glueball. TH models moreover contain all the structure required to incorporate dark matter (DM) candidates, in the form of either a twin weakly charged WIMP [6,7] -hereafter, a 'T-WIMP' -or twin baryonic asymmetric dark matter [8,9]. In this work we focus on TH models featuring T-WIMP candidates.
Although they require no additional light states with SM charges, twin Higgs models may nevertheless produce striking signatures in a variety of direct and indirect searches.
The effectiveness of the various probes is mostly determined by the lifetime, decay modes and masses of the LTH and LTP, and as such it is useful to sketch out the space of signatures in terms of the properties of these two particles. In detail: (i) Relativistic LTPs at either the BBN or CMB epochs lead to strong tensions with current bounds on the effective SM neutrino degrees of freedom ∆N eff ; (ii) Twin states may be produced at the LHC through an exotic decay of the SM-like Higgs or its heavier partner, or through the decay of some of the TeV-scale, colored states that could be accessible in certain UV completions. In some TH models, this can lead to interesting LHC signatures from the decay of the LTH [5, [10][11][12]; (iii) A LTH that is metastable on BBN timescales may lead to either overclosure, arising from other (meta)stable hadrons with masses nearby to the LTH, or post-BBN matter-dominated eras that disrupt the standard BBN paradigm; (iv) Annihilation of T-WIMP candidates into hard b quarks, charm quarks or τ 's may create tensions with evolving astrophysical bounds on antiproton or positron production in the galactic central stellar cluster (see e.g. [13][14][15][16][17][18]).
Finally, apart from these LTH and LTP signatures, one may also search for deviations of Higgs couplings. These are in principle the most robust probes of twin Higgs models, but their reach at the LHC is rather limited [19].
Several TH models have been previously proposed, that may be broadly characterized by these signatures. The first instance is the 'mirror twin Higgs' [1], which includes a full copy of the SM in the twin sector. The LTPs are therefore twin neutrinos and twin photons, while the LTH is a twin pion that rapidly decays to twin leptons and twin photons well before BBN. This model is susceptible to strong BBN and CMB bounds on ∆N eff . To evade these, one either requires an asymmetric reheating between both sectors [20][21][22] or require large entropy production during the QCD phase transition [1,23,24].
Two other recent TH proposals include the 'fraternal twin Higgs' [5] and the 'vectorlike twin Higgs' [25]. In these constructions, the low energy spectrum of the twin hadronic sector is given by a zero or one light flavor twin QCD theory. The LTH is then either a glueball or an onium state whose mass is determined by the twin confinement scale, Λ qcd , that is itself constrained to be nearby the SM confinement scale, Λ qcd . In the vector-like twin Higgs model, the LTH and the LTP are the same particle, but in the fraternal twin Higgs model the LTP can be either a twin neutrino or the LTH itself. The decays of the LTH may generate (possibly displaced) collider signatures, while a twin neutrino LTP contributes ∆N eff 0.075 [6,7], that is potentially detectable in the future. These models can further admit a twin tau T-WIMP that annihilates into twin hadrons, some of which in turn decay to bb or ττ pairs (see e.g. [7]), potentially detectable in (future) astrophysical data.
In Table I we schematically summarize the signature space of these twin Higgs models according to the signatures and potential features (i)-(iv) above. Laid out in this fashion, it becomes clear that there is a fourth category of TH models, which is so far largely unexplored in the Literature. Namely, the case where the LTP is heavy enough to trivially avoid ∆N eff constraints, while the LTH is at the same time always light enough, and therefore metastable enough, to be invisible at the LHC. While certain parts of the parameter space of the fraternal  and vector-like twins realize this scenario, we consider here models in which these features are robust predictions everywhere in parameter space. A well-motivated and representative example is a twin sector that contains a mirror copy of the SM hadrons -i.e. all three twin quark generations with three light quark flavors -but no dark radiation -i.e. no light leptons or photons. A T-WIMP may be incorporated as an additional Dirac or Majorana state. We call this scenario the 'hadrosymmetric twin Higgs'.
In this scenario, the existence of multiple light twin quark flavors produces a twin pion LTH, a pseudo-Nambu-Goldstone boson (pNGB). Unlike in zero and one light flavor twin sectors, here the pion mass may be much lighter than the twin confinement scale. The absence of dark radiation ensures the twin pion is also the LTP and avoids, in the most straightforward manner, tension with ∆N eff bounds, provided the twin pion is heavier than T BBN ∼ MeV. Without dark radiation decay modes, the twin pion must promptly decay by BBN to SM degrees of freedom to avoid overclosure or a matter-dominated era thereafter by twin pNGBs. On collider timescales the twin pion is then stable (cf. e.g. Refs [26,27]) and appears as missing energy produced with only a small cross-section: From the viewpoint of collider and ∆N eff searches, the hadrosymmetric scenario then represents a 'least detectable' scenario compared to the other models in Table I (This situation is similar to the case with one light flavor, where the LTH is 0 −+ meson [7,8].) There are a limited number of choices for such a portal, especially if one wants to retain a minimal twin-Higgs sector and not introduce additional SM-charged particles. There are two phenomenologically different regimes for this portal, corresponding to whether the twin pion mass is above or below three times the SM pion mass, 3m π : Since the twin pNGBs can be much lighter than the twin confinement scale, the twin pion mass is a free parameter of the theory. In the case that the twin pion mass is above 3m π , we find the portal does not require a UV completion below the twin Higgs cutoff scale. However, if the twin pion mass is below 3m π , we show, for a representative mediation model, that in most of the viable parameter space the BBN lifetime bound implies that the mediator itself should be accessible either at the LHC or at the intensity frontier.
Second, the strongly-coupled twin sector generically acts like a neutral natural hidden valley [28,29], and in the spirit of Ref. [30], DM annihilations into twin quarks can therefore produce 'dark showers' of twin hadrons, that then subsequently decay to the SM sector. In the case that the twin pions are lighter than 3m π , a leptophobic mediator ensures that twin pion decays are predominantly to diphotons. If instead they are heavier than 3m π , the twin pions decay mainly to three SM pions and hence to a higher multiplicity of softer photons or leptons. The indirect detection signatures for dark showers are then mainly encoded into photons, as in Ref. [30], evading present or future astrophysical cosmic-ray bounds on the associated emission of antiprotons and positrons.
This dark showering to photons may produce a visible astrophysical signal. We present detailed numerical simulations of the photon spectra produced by Dirac and Majorana T-WIMP annihilation in the galactic center, compared against γ-ray data from Fermi -LAT.
We find that for a Dirac T-WIMP, there is a viable region of parameter space in which these dark showers may reproduce the putative galactic center γ-ray excess (GCE), claimed to be seen in Fermi -LAT data from the central regions of the Milky Way galaxy [31][32][33][34][35][36][37][38][39] and dwarf spheroidal galaxies [40,41], as well as simultaneously generate an appropriate amount of DM and exhibit minimal fine-tuning. (See also Refs. [39,[42][43][44][45][46][47][48] which discuss alternative astrophysical explanations for the GCE, or dispute observations of an excess in dwarf spheroidal galaxies [49].) Alternatively, treating observed photon fluxes as a conservative upper bound, there are large regions of parameter space are not excluded by current galactic center γ-ray data.
This paper is organized as follows. In Sec. II we present the hadrosymmetric twin Higgs model and benchmark dark matter models, followed by an examination of all applicable astrophysical, cosmological and collider bounds in Sec. III. In Sec. IV we consider effective field theory analyses for portals mediating the twin pion decay, as well as a sample UV completion and corresponding constraints on its mediator. Simulations of dark showering DM annihilations and corresponding photon spectra are presented in Sec. V, along with both fits to GCE spectra and constraints from γ-ray fluxes.

A. Twin Higgs framework
In the twin Higgs framework the SM Higgs is realized as the pNGB of a spontaneously broken global symmetry. The key difference with traditional pNGB Higgs frameworks is that the symmetry in question is accidental rather than exact. This has the advantage that the top partner may be fully neutral under the SM gauge groups, which removes conventional collider constraints on the top partner. On the other hand, this also implies that the accidental symmetry is generally not radiatively stable beyond 5 to 10 TeV, at which scale a UV completion is then required. The known options for UV completions are supersymmetry [24,50,51], compositeness [52][53][54], orbifolds [3,4] and holographic setups [55,56].
where y t (y t ) is the top (twin top) yukawa, and g 3 (g 3 ) is the (twin) color gauge coupling.
There are no restrictions on the twin hypercharge coupling, g 1 , nor on the remaining twin yukawa couplings, other than that they cannot be O(1). The complete absence of the first two generations, or a gauged twin hypercharge, do not reintroduce an appreciable amount of fine-tuning. In this paper we consider the 'hadrosymmetric' twin Higgs scenario, for which the twin sector contains twin copies for all three generations of quarks, but the twin lepton sector is absent. We also assume that twin hypercharge, if present, is at most only a global symmetry, but is not gauged. In the absence of gauged twin hypercharge, the heavy SU (2) L gauge bosons are degenerate, m Z ,W = g 2 f /2, so that their masses scale with a factor of f /v with respect to their SM counterparts. Naïvely, all twin quarks are similarly typically f /v times heavier than their SM counterparts. However since all twin yukawas except y t are almost unconstrained, this assumption may be relaxed for all quarks except for the twin top. For most of the discussion in this paper, we nevertheless assume the twin yukawas are set by the SM-twin Z 2 such that y = y, and the twin quark masses are then for all q .

B. Twin hadron spectrum
For the O(10%) deviation of twin QCD coupling, g 3 , permitted in eq. (2.1), computing the two-loop renormalization group evolution from the UV cutoff down, one finds a corresponding O(1) deviation in the twin confinement scale compared to the SM, viz.
With two or three light flavors, the dependence of λ on the mass spectrum of the twin fermions is rather mild. We shall therefore treat the confinement scale ratio λ as a free parameter of the twin theory, within the range indicated by eq. (2.3).
We characterize twin quarks as light or heavy depending on whether their mass is small  Fig. 1 we show this approximate region.
As in the SM sector, the lightest pNGBs, and thus the lightest twin hadrons, are twin pions, followed by heavier twin kaons and η's. The twin pions embed into a twin isospin triplet (π ± , π 0 ). Unlike in the SM sector, in the absence of a gauged twin hypercharge, twin isospin breaking arises only from the splitting of the yukawas. In particular, twin pion mass splitting effects arise only at second order in isospin breaking, such that in which δ SM ∼ 1% is the isospin breaking in the SM sector induced by splitting of the up-down yukawas. In most of the parameter space, we therefore expect the mass splitting of the π 0 and π ± to be much smaller than for the SM pions, but with the same sign.
In the class of twin sector scenarios we consider in this paper, the twin pions π ± are charged under a (residual) global U (1) and therefore stable. In the absence of a light twin photon or other twin-SM mediator, the π 0 can only decay via the Higgs portal to SM final states. This decay requires both parity and twin isospin violation and is therefore heavily suppressed. We discuss this and twin-SM mediator models for decay of the twin pion in further detail in Sec. IV.
Apart from pions, the remainder of the twin light hadron spectrum comprises: higher spin states, such as the ρ or ω ; a twin η ; glueballs such as the f 0 ; and baryons, such as the twin proton p and neutron n . Similar expectations apply to the expected twin protonneutron mass splitting, as for the pion mass splitting. The twin proton is stabilized by baryon number, while the heavier twin neutron is stable as well due to the absence of light twin leptons. We assume all heavier hadrons strongly or weakly decay to pions, n and p , as in the SM sector, and that there is no twin baryon asymmetry.

C. Dark matter
We consider two different twin sector benchmark scenarios that produce respectively Majorana or Dirac DM, coupled to the twin sector by twin electroweak interactions. Along with three generations of twin quarks -electroweak doublets Q i and singlets u c i and d c i , for i = 1, 2, 3 -we restore a single lepton-like generation -an electroweak doublet L and a singlet E -and an anomalous global twin hypercharge symmetry. This single electroweak doublet is the minimal content required to cancel the SU (2) L Witten anomaly. One then obtains an effective theory that may be UV completed without anomalous global symmetries.

A. Relic density
The thermally averaged 2 → 2 annihilation cross-section for fermions of mass m DM at Here S is a symmetry factor -S = 2 for Majorana χ and S = 1 for Dirac ψ -and Σ(s) is the spin-summed square amplitude for annihilation to two (twin) quarks, integrated over the final state phase space. Including both annihilation to twin quarks through the twin Z as well as annihilation to SM bb via the Higgs portal, one finds where m Z = m Z f /v is the twin Z mass and we treat the twin quarks as massless, so that annihilation through the twin Z proceeds to N f = 5 twin quark flavors. The latter assumption, which greatly simplifies the square amplitude expressions, anticipates the result that m χ,ψ m b for Ω χ,ψ = Ω DM ; the regime that m χ,ψ < m b cannot achieve a sufficient DM relic abundance. Corrections to σv on the Ω χ,ψ = Ω DM contour arising from the b mass arise at the O(%) level, and can be safely neglected. From the common s − 4m 2 χ factors, observe that the Majorana cross-section is p-wave suppressed, as expected from Fermi statistics. Comparing to the fraternal twin Higgs, the annihilation cross section is somewhat enhanced here because of the larger number (N f = 5) of light twin quarks available as final states. The required DM mass will therefore be somewhat smaller than was obtained for the fraternal twin WIMPs [6,7].
To a good approximation, one finds that the thermal relic abundances in which the Higgs resonance contributions are not displayed; x f = m DM /T f , the freeze-out temperature; the functions and the parameter We take the form factor g HN 1.2 × 10 −3 for both proton and neutron from lattice studies [62,63]. In Fig. 3  below the current LUX bound [64]. Anticipated future sensitivities [65] may, however, have sufficient reach to probe this scattering.

E. BBN bounds
The twin and SM sectors typically decouple at the GeV temperature epoch, before confinement of the strong sectors. Without a twin baryon asymmetry, and without light twin leptons or radiation, at late times the entropy and energy density of the twin sector is deposited dominantly into the lightest hadron, here the π 0 . In order to avoid a π 0 matterdominated phase at or beyond the BBN epoch, or overclosure by π ± , we therefore require the twin pion to decay into SM degrees of freedom before BBN (see similarly [8]), i.e.
One may be concerned that the stable π ± may still freeze out from the hadronic plasma with a significant relic abundance, since the isospin-breaking induced mass splitting with the π 0 is small. We therefore now verify that the π ± relic abundance is negligible.
mann equations for π ± freeze out are Here we assume the yield Y + = Y − , and define the mean square pion massm 2 π ≡ (m 2 π + + m 2 π 0 )/2 and x ≡m π /T . The pion-pion scattering matrix element (s − m 2 π )/f 2 π [66], whence one may show the forward and inverse thermally averaged cross-sections Fig. 4 we show the thermal history of the charged twin pion relic abundances for a range of mass splittings ε. Even in the ε 1 regime, the relic abundance for the π ± is Ω π ± ∼ 10 −5 Ω DM , so there is no significant twin pion thermal relic. During the recombination epoch, charged twin pions with this small of an abundance annihilate into SM photons with thermal cross section σv ∼ 10 −32 cm 3 s −1 . This is well below the current upper bound ∼ 10 −29 cm −3 s −1 for m π 100 MeV, arising from measurements of cosmic microwave background (CMB) anisotropies [60,61].
Note moreover that for a twin pion of mass m π ∼ 100 MeV, number changing processes such as 4π → 2π remain efficient until m π /T ∼ 10 (see e.g. [67]). If the twin sector is kinetically decoupled from the SM, then such interactions may cause the twin sector temperature to rise exponentially compared to the SM sector. This keeps the π ± in equilibrium with π 0 until later times, further suppressing the charged pion relic abundance.

IV. TWIN PION DECAY
In the absence of twin hypercharge interactions, the twin pion can nominally only decay by mixing with a 0 ++ isospin singlet state, which can subsequently decay through the Higgs portal to SM degrees of freedom. Such mixing is however heavily suppressed by both parity and the small twin isospin violating Yuwaka coupling, y u − y d . For instance, the twinisospin and parity-violating π 0 G µν G µν /f π operator enters only at least at two loop order [68], since the only source of parity violation is the twin electroweak interactions, assuming a negligible θ qcd . The twin pion decay amplitude is thus heavily suppressed, leading to lifetimes well in excess of τ BBN . In order to allow the π 0 to decay before BBN, it is therefore necessary to introduce an additional portal between the twin and SM sectors, In this section, we first perform an effective operator analysis for such interactions. For m π 0 < 3m π 0 , this analysis reveals the need for a UV completion of the additional twin-SM portal below the TeV scale. We subsequently present a sample UV completion and analyze the corresponding constraints on its parameter space, both from direct searches at the LHC as well as from precision flavor measurements. Alternatively, in the regime that 3m π 0 < m π 0 < 2m K , hadronic decays of the twin pion to three pion final states -i.e. π 0 → 3π -predominantly produce soft photons, and are sufficiently fast that no extra twin-SM mediator is required below the twin Higgs cutoff, Λ ∼ 5 TeV. Hence in this regime the hadrosymmetric twin Higgs framework does not require any near-term detectable phenomenology in order for the twin pions to decay before BBN. In the regime m π 0 > 2m K , the available SM hadronic decay channels of the twin pion become large in number. This case would also require a rather large f /v 10 from eq. (2.6), and we therefore do not consider it in this work.

A. Effective field theory analysis
The twin and SM confinement scales are readily probed at current collider experiments, so we therefore must consider a (twin) quark effective theory for the twin pion decay. The twin pion decay must be mediated by a twin isospin-violating current, to avoid suppression of the decay rate by the small y u − y d coupling. Under the approximate Z 2 symmetry between the SM and twin quark sectors, this current should similarly couple to SM quarks rather than leptons, and be SM isospin-violating, too. In this framework, then, the π 0 decay is maximal in the case that it mixes with the lightest SM pseudoscalar in a non-trivial isospin representation, that is, the π 0 .
Gauge kinetic mixing portals with SM hypercharge cannot induce π 0 tree-level diphoton or dilepton decays, because they do not mix longitudinal modes of heavy vector bosons, and therefore generate an insufficiently fast twin pion decay. The lowest dimensional viable portals are then dimension-six current-current interactions. We consider vectorial or axial neutral currents, and further require them to be flavor-conserving to avoid precision flavor constraints. Portals involving the V −A current Q † σ µ Q do not generate the required isospin breaking in the twin sector. We therefore focus only on V + A interactions, which couple to both SM and twin right-handed quarks.
Defining the current for right-handed q = u i ,d i ,u i or d i , pion-twin-pion mixing is generated by the effective where Λ P is the scale of the twin-SM portal. This operator is parity violating but CP conserving, and therefore permits a J P C = 0 −+ state to mix with either 0 −+ or 0 +− states.
Before proceeding, it should be mentioned that there are several alternative portals, other than a current-current operator. First, one might generate the π 0 decay with a spin-0 mediator coupling to Q u c . Such a mediator must belong to a twin electroweak doublet.
This could be implemented by extending both the Higgs and twin Higgs sectors, for example in a twin two Higgs doublet model. However, since this additional Higgs doublet must couple to the quark sector, the Glashow-Weinberg condition is violated and one generically expects large deviations in precision flavor experiments. While this may be interesting direction for further work, we do not pursue this option further here. Second, one could instead consider explicitly breaking the SM-twin Z 2 exchange symmetry, by allowing a π 0 decay channel to leptons, mediated by the effective operator where = µ or e for the π 0 masses under consideration herein.

B. Decay rate and branching ratio estimates
The amplitude for π 0 → SM decay, as generated by the operator (4.2), receives its dominant contribution from the off-shell π 0 channel, i.e π 0 → π 0 * → SM. The next lightest isospin triplet pseudoscalar is the much heavier π(1300), and we neglect Yukawa-suppressed SM isospin violating effects. To estimate the decay rates, observe that the mixing amplitude in which we have applied the scalings in eqs. (2.4) and (2.6). Thus the decay rate and branching ratios for the π 0 are immediately informed by estimating the corresponding partial widths of an off-shell π 0 * , with mass m π 0 .
It is therefore not mediated by the CP-conserving operator (4.2), and is negligible. Finally, for m π 0 > 3m π 0 , the twin pion can decay to three SM pions and further exclusive hadronic modes. These strong decays are not suppressed by isospin violation, and are therefore expected to be much faster than the electromagnetic decay modes.
To estimate the branching ratios of the twin pion with eq. (4.5), we adapt data from the SM π 0 and η -the neutral pseudoscalars closest in mass to the m π 0 regime of interest -rescaling the partial widths according to their explicit mass dependence. To adapt the η data, we take care to also rescale by Clebsch-Gordan coefficients associated to decays of an isospin triplet rather than a singlet, as well as including isospin violating effects and η-η mixing. (We emphasize that the η here is the SM hadron, not a twin state.) One finds, in a particular basis of isospin invariants, that where φ is the η-η mixing angle, δ SM ∼ 1% is the SM isospin breaking scale, and r 1,2,3 are unknown ratios of hadronic matrix elements. Estimates provide φ 40 • (see e.g. [69]). (For the ideal mixing case that the η is a pure ss state, φ = tan −1 √ 2 55 • .) Taking a naïve average over O(1) values for r 1,2,3 , for 2m π < m π 0 < 3m π one estimates In eq. (4.7) we assume that the η → π + π − γ process is vector-meson dominated, so that it can be thought of as η → (ρ 0 * → π + π − )γ, and hence its rate naïvely scales as m 7 η . We do not include the phase space effects that smoothly turn off the π 0 → π + π − γ rate nearby to the m π 0 = 2m π threshold. For m π 0 > 3m π , similarly one estimates where the 3π state is either π + π − π 0 or 3π 0 . The η → 3π decay naïvely scales as m 5 η /f 4 η , whence the mass scaling dependence in eq. (4.8). Note that relative Clebsch-Gordan coefficients, combinatorics and phase space symmetry factors imply that Γ[π 0 * /η * → 3/2. Thus in the π 0 → 3π decay mode, one expects the photon versus lepton production ratio to be approximately 11 : 4. Therefore, in the regime 3m π 0 < m π 0 < 2m K one expects twin pion decays to predominantly produce a high multiplicity of comparatively soft photons.
The corresponding estimated twin pion branching ratios as a function of m π 0 are shown in the left-hand panel of Fig. 5, for f π = f π . We see there that for m π 0 < 2m π , the diphoton rate dominates as expected. For 2m π < m π 0 < 3m π the diphoton mode continues to dominate the π + π − γ decay mode, while for m π 0 > 3m π , the purely hadronic 3π modes π → γγ π → π + π − γ π → 3π 0 π → π + π − π 0 0.1 dominate overwhelmingly. In the righthand panel of Fig. 5 we show the maximum effective scale Λ P under which the twin pion decays before BBN -i.e. τ π 0 1 s -in these different kinematic regimes. We see in Fig. 5 that to satisfy the BBN bound for m π 0 < 3m π , the portal (4.2) generically requires a UV completion nearby the TeV scale, while for m π 0 > 3m π , the BBN bound is satisfied generically for Λ P well above the scale of twin Higgs effective theory, Λ ∼ 5 to 10 TeV. Since the LHC probes TeV energies, an effective operator approach is therefore insufficient to study possible collider constraints for m π 0 < 3m π . We shall therefore present a sample UV completion for this regime in the next section, as well as its corresponding experimental constraints.
Finally, if the leptonic portal (4.3) is available, then the twin pion may also decay rapidly to leptons, subject to the usual chiral suppression. In particular, whence one may compute the decay rate directly, viz. (4.10) Hence if kinematically accessible, the muon decay channel dominates, followed by a much weaker decay to electrons. In the righthand panel of Fig. 5 we show also the maximum effective scale Λ P required to ensure τ π 0 1 s via the leptonic portal (4.3). The muon decay channel is fast enough by itself that the twin pion lifetime BBN bound can be satisfied with Λ P 10 TeV. Therefore, in this scenario no UV completion is required below the twin Higgs cutoff. (Alternatively, one could also modify the model to include gauged twin hypercharge, such that the twin pion can decay to two twin photons [8]. The twin photons then subsequently decay into SM leptons through the kinetic mixing portal.) However, this leptonic decay channel may introduce possible tensions with astrophysical bounds on high multiplicity muon production by DM annihilations [13,14]. A full analysis of these astrophysical bounds is beyond the scope of this work, so we shall therefore not consider this channel further.

C. A sample UV completion
A straightforward UV completion of the operator (4.2) can be achieved by charging the right-handed quarks under an additional broken U (1) X gauge interaction, such that u (u ) and d (d ) have opposite charges. 2 That is, an interaction of the form where X is the massive U (1) X gauge boson, and g X is the gauge coupling. Note also that charging the right-handed quarks in this manner introduces mixed U (1) X -hypercharge instanton anomalies. However these can be cancelled with a small set of TeV-scale hypercharged anomalons.
To accommodate the BBN bound on the twin pion lifetime for m π 0 < 3m π , from Fig. 5 one requires m X /g X = Λ P 2 TeV . (4.12) We consider the range 10 MeV < m X < 2 TeV and allow g X to vary to accommodate this BBN constraint. 2 Here we choose to adhere to the spirit of neutral naturalness by not introducing additional, TeV-scale colored states. There exist, however, UV completions of the operator (4.2) if we lift this requirement.
For example, a t-channel scalar interaction of the form in which the Yukawas ξ may be dynamically aligned with the SM and twin Yukawas (see e.g. [70]). The scalars φ must carry both twin and SM colors and therefore this UV completion suffers strong constraints from collider searches. Options of this type were also explored in a different, emerging-jet context in [26].
The U (1) X charge assignments imply that the twin and SM Yukawas must be modified to the formỹ where φ is a scalar charged under the U (1) X , and the interaction scales Λ u,d φ may be different for up and down-type couplings, but are at least as high as the twin Higgs cut-off, Λ. The SM Yukawas are then generated by y u,d =ỹ u,d φ /Λ u,d φ . The vev φ contributes to the mass of the X, but may not be the only contribution in a complete model, and we therefore keep the mass of the X, m X , as an independent parameter.
If the X couples universally to all three quark generations, the requirement of an O(1) TeV, there will be a region of parameter space with relatively large g X and low m X where some additional fine-tuning is needed to ensure that the X remains light compared to g X φ . But, as we will show in the following section, the tuned region of parameter space for m X > 1 GeV happens to be outside the experimentally allowed regions, and therefore this tuning is not of great concern.
An alternative is to imagine that the U (1) X does not couple to the third quark generation, One might have also considered charging the Higgs under the U (1) X , thereby eliminating the need for the spurion φ in the quark sector. (Since the leptons are uncharged under U (1) X , a spurion would still be needed for the lepton Yukawas.) This approach is, however, plagued with problems, as the Higgs vev now induces a tree-level mixing between the X and both the SM Z and twin Z bosons. This induces strong tensions with both electroweak precision tests and DM direct detection experiments. In particular for the case where the dark matter is a Dirac fermion, ψ, the mixing induces the four-fermi operator ψ † γ µ ψ q † γ µ q/Λ 2 P . This is ruled out by current direct detection experiments, given the BBN bound (4.12). We therefore do not consider this approach in this work.

D. Model constraints
At high m X -above several hundred GeV -the only relevant collider constraints come from LHC searches, while for low and intermediate m X -below and above 1 GeV respectively -a variety of experiments at the intensity frontier play a role as well. In what follows, we separately consider the mass ranges m X > 300 GeV, 300 GeV > m X > 1 GeV and m X < 1 GeV, at the following benchmark points: benchmark 1: m π 0 = 250 MeV, f π = f π , benchmark 2: m π 0 = 350 MeV, f π = 2f π , (4.14) where the latter benchmark anticipates the best fit value for the galactic center γ-ray excess, discussed in Sec. V. m X > 300 GeV: Rescaling the dijet and monojet LHC bounds [72], one finds that m X > 500 GeV is presently already directly excluded. For lower masses and lower couplings these searches lose sensitivity, but the combination of the monojet bound with the BBN constraint on the twin pion lifetime still excludes m X 300 GeV.
1 GeV < m X < 300 GeV: In this regime LHC monojet searches as well as a variety of dark photon searches at B factories apply, but the BBN bound is the only bound which depends on the choice of {m π 0 , f π }. In Fig. 6   In the gray shaded region for the universal case, at least 10% fine-tuning is required to keep the X light.
where we take the branching ratio of X to twin quarks to be 50%. The most important cut for this search on the missing transverse momentum / E T > 350 GeV (see Ref. [73] for more details). With the same background analysis as in Ref. [78], we obtain a 95% CL upper bound on g X , which has only a very mild dependence on m X . For the 14 TeV projection, shown by the red-dashed curve, we keep all the cuts in the 8 TeV search except / E T > 550 GeV. From the projection of a data-driven analysis [78], we obtain this bound with 300 fb −1 of data and an estimated 3.4% systematic uncertainty. Unless there are significant improvements in the reduction of uncertainties for the signal acceptance and selection efficiency, the systematic uncertainty on the 14 TeV search is not expected to be much lower than the CMS study.
Hence the improvement of this bound is very limited compared to the 8 TeV search. The dijet constraints in this region of parameter space come from the UA2 experiment [79] but are always much weaker than the monojet constraint. We therefore do not include them in  (4.15) assuming = 0 at the twin Higgs cutoff scale, Λ = 5 TeV. Here g 1 is the SM hypercharge coupling, and Y f and q f are the hypercharge and U (1) X charge of the fermion f respectively. 4 In the universal X coupling case, in which X couples to the third generation, this can be approximated by while the analogous expression for the horizontal case is As a result of this kinetic mixing, the SM leptons are millicharged under the X. Note that because of the tree-level couplings of X to SM and twin quarks, the X branching ratios to jets or missing energy nevertheless dominate leptonic branching ratios when m X > 2m π .
The consequent electroweak precision (EWPT) bounds in Fig. 6 were rescaled from those obtained in Ref. [80]. The dilepton constraints below m Z are rescaled from those in Ref. [81], where we use Madgraph 5 [82] and Feynrules 2.0 [77] to compute the cross sections. The projected bound, shown by the dashed orange line, assumes 3 ab −1 of 14 TeV data. The dilepton bound above m Z is obtained from a similar recast of the CMS search for A 0 → µ + µ − [83]. Since the X signal does not contain any b-tags, we only make use of 'category 3' events in Ref. [83], which corresponds to gluon-fusion production. We further implement a correction factor of 0.7 to account for the difference in acceptance between this search and the X signal. Given that the branching ratio of the X to leptons is very small, there is no meaningful bound from h → ZX → 4 [80].
There are moreover several applicable B factory searches, that probe this effective kinetic mixing. Firstly, we rescale the result for the search for e + e − → Xγ that was carried out in Ref. [84]. The 95% CL excluded region is shown in brown in Fig. 6. The dashed brown line indicates the expected exclusion reach from Belle II, also rescaled from Ref. [84], using a projected 50 ab −1 of data at the Υ(4S) resonance. Secondly, we adapt the branching ratio B(Υ(1S) → invisible) < 3×10 −4 at 90% CL from BaBar [85]. In the universal coupling case, this can be interpreted as a bound on the invisible decay Υ(1S) → X * → q q decay, since for m X ≥ 2m π and g X 0.1 the twin hadrons produced in the X decay all decay outside the detector. We extract a bound on g X from Refs. [86,87] by comparing the invisible branching ratio with the well-measured branching ratio of Υ(1S) → µ + µ − , where we assumed that the unknown hadronic matrix elements for each process are approximately the same. The result is shown in the purple-shaded region in Fig. 6. The projected sensitivity of Belle II is indicated by the dashed purple line, and is derived by assuming the branching ratio constraint (4.18) will improve to the 4 × 10 −4 level, as discussed in Ref. [88]. Finally, for the universal case it is plausible that a slightly stronger bound can be obtained by recasting the BaBar search for Υ(3S) → γ+pseudoscalar [89]. However, in this case a careful treatment of the hadronic matrix element is necessary, which is beyond the scope of this paper. (See for instance Ref. [86] for an analysis in terms SM-DM effective operators.) We see in Fig. 6 that the projected reach of Belle II dark photon searches and dilepton searches should probe nearly all of the presently allowed parameter space. m X < 1 GeV: In this regime a different set of constraints becomes relevant. First, the EFT approximation leading to the BBN bound (4.12) is no longer valid, and the X must be included as a dynamical degree of freedom. This does not significantly affect the bound, except for m π 0 nearby to the X resonance or if m X < m π 0 /2. In this latter case the twin pion can directly decay to two on-shell X bosons through the U (1) X chiral anomaly, rather than a diphoton final state, and the BBN bound no longer depends on the X mass.
Second, since α s is large in this regime, one expects O(1) incalculable corrections to the mixing (4.15), which increases the theoretical uncertainty on the bounds that rely on kinetic mixing. Here we continue to use eq. (4.15) for definitiveness, while keeping these O(1) uncertainties in mind. In Fig. 7 we show all applicable constraints in this regime for the universal model. The horizontal model is nearly identical. The BaBar and monojet bounds, indicated by the brown and red shaded region respectively, are sensitive to various thresholds involving SM and twin pion masses. Specifically, once m X < 2m π , the X can no longer decay to the twin sector. As a result both the monojet bound and the BaBar bound on γ+MET disappear. In this regime the most powerful BaBar bound comes from X → + − decay [90], which begins to dominate once the decay to SM pions is forbidden, i.e. for m X < 2m π .
Finally, for very a light X there are constraints from π → Xγ decay at the NA48/2 experiment [91], as well as from beam dump experiments [92]. For completeness we also include projected bounds from APEX [93] and HPS [94], in addition to a projected bound from D * 0 → D 0 X at LHCb [95]. Since the X couples directly to the quark sector, the rate in NA48/2 and LHCb receives a contribution in addition to that from the kinetic mixing.
Since this contribution comes with an incalculable hadronic matrix element, we simply take it to be zero here, which is very conservative.
In summary, the constraints are qualitatively similar for the universal and horizontal cases, and that much of the remaining available parameter space will be probed in upcoming experiments.
Since the π 0 s decay mostly to photons -diphoton final states for m π 0 < 3m π , and higher multiplicity soft photons for 3m π < m π 0 < 2m K -the dark shower of light twin hadrons produced by the DM annihilation can result in an astrophysically detectable signal. Compared to stable (or long-lived) charged decay products, the production spectra and distribution of photons do not suffer from large corrections arising from propagation, and therefore provide robust probes or constraints on the underlying DM distribution and annihilation process.
Over the last decade, much attention has been focused on a galactic center γ-ray excess (GCE), claimed to be seen in Fermi -LAT data from the central regions of the Milky Way galaxy [31][32][33][34][35][36][37][38][39] that may be a signal of DM annihilations into a variety of SM final states.
Corresponding excesses in dwarf spheroidal galaxies are also claimed to be seen [40,41], though this is disputed by the Fermi collaboration [49]. The observed GCE flux corresponds to an underlying DM annihilation rate comparable to a thermal relic cross section, and exhibits a spherically symmetric morphology consistent with a NFW-type DM profile.
Although the DM explanation of the GCE has been called into question [42][43][44][45][46][47][48], and should be approached with caution, models with DM annihilations at rates comparable to thermal relic cross-sections can generically produce photon fluxes comparable to current experimental sensitivities of the Fermi -LAT. Hence, any thermal DM model that predicts γ-ray signals from annihilation needs to be consistent with the observed γ-ray spectra, no matter if treated as a signal or constraint.
Cosmic rays -antiprotons and positions -are typically produced alongside γ-ray spectra in many frameworks designed to address the GCE, possibly resulting in mild tensions with Pamela and AMS-02 data (see e.g. [13][14][15][16][17][18]). However, DM annihilation into cascade processes provides a generic way to soften the spectra of cosmic rays, and hence loosens the corresponding astrophysical bounds [96]. Most previous analyses have focused on cascades with fixed multiplicity, but another option is to consider showering in a strongly-coupled dark sector [30]. This is precisely the situation arising from T-WIMP annihilations into showers of light twin hadrons within the hadrosymmetric twin Higgs framework. In this section we adapt our prior treatment in Ref. [30] of this so-called dark showering to simulate the γ-ray spectra produced by such DM annihilations.
The spectra are generated by Pythia8 [75], with the QCD parton shower and hadroniza-tion model modified appropriately for the twin quark sector. In detail, in order to be able to correctly approximate the twin shower spectrum, the running of Pythia8 has to be suitably generalized to allow for the variation of twin quark masses relative to the twin confinement scale. This involves two considerations: First, the mass thresholds for the running of strong coupling and the parton shower cutoff are adjusted relative to the quark masses to make sure that hadronization into the heaviest pNGBs of the theory -the analogs of the K and η mesons -is still kinematically allowed. Second, the spectra of hadrons has to be appro-  However, these heavier states have a less significant effect on the spectrum of the dark shower.
Computation of astrophysical rates is done using PPPC 4 DM ID [97]. With this setup, we proceed to assess the compatibility of our DM annihilation signals with the observed Milky Way galactic center γ-ray spectrum, whether treated as an actual spectrum to be fit, or as a possible background providing a constraint.

A. GCE fits
In order to understand the effect of the twin parton shower on the γ-ray spectrum, we Goodness-of-fit (χ 2 /dof) contours for fits of photon spectra produced by T-WIMP annihilations to the GCE signal are shown in Fig. 8, both for m π -Λ qcd and f /v-λ parameter spaces, while the spectrum at the point of best fit is shown in Fig. 9, and compared with the best fit for direct annihilation into bb pairs. Because of correlations between the energy bins [38], the best fit curves differ from those that would be expected under the assumption of purely uncorrelated signal bins. Namely, while the overall normalization is dominated by the lowest energy bins, the shape of the spectrum is predominately determined by the peak and the beginning of the tail region. This results in a fit that appears to anomalously undershoot the highest-energy bins.
Since the uniformity of the scaling relations (2.4) and (2.6) should only be taken as approximate, we have checked to see the effect of varying the masses of the excited multiplets by an additional 10% around their nominal values. As twin quark masses relative to the confinement scale increase, this variation introduces larger and larger deviations in the spectrum. In Fig. 8, we shade in red the region where these deviations change the χ 2 /dof goodness-of-fit metric by more than unity, which we take as a sign that the shower is sensitive to the non-perturbative aspects of the hadron spectrum at a level that makes spectrum  Red shading indicates regions of large sensitivity to the details of the twin mass hadron spectrum (see text for more details). Spectra provided by Ref. [38].
The fitted DM masses range only from 56 to 59 GeV over the entire displayed parameter space for both scenarios, with the Dirac and Majorana DM best-fit masses both being 57 GeV. In both cases, the best GCE fit regions are partially excluded by Higgs invisible width and CMB reionization bounds. For the fits shown in Fig. 8, we can see that the best fit regions for the spectra correspond to a twin confinement scale slightly higher than that of the SM, and the twin quark mass to confinement ratio m q /Λ qcd = (f /v)/λ(m q /Λ qcd ) is slightly higher than in the SM sector, too. Moreover, the best fit region lies in the three In Fig. 10 we show the GCE goodness-of-fit contours together with contours of constant DM annihilation cross-section as extracted from the GCE spectral fits, rather than from the relic abundance theory contours (3.3). In the Dirac case, the annihilation cross-section contour corresponding to the observed DM relic density, σv 3 × 10 −26 cm 3 s −1 , intersects the best fit GCE region, which is by itself a non-trivial agreement. We see further that this annihilation cross-section contour moreover intersects the Ω ψ h 2 = 0.12 theory contour in the best fit region: A remarkable concordance between predictions of the hadrosymmetric twin Higgs framework and the observed GCE spectral features.
Besides GCE signals, DM capture and subsequent annihilation in the Sun might also provide detectable γ-ray spectra. In the Dirac scenario, ∼ 60 GeV DM particles have a nuclear scattering cross-section σ N 10 −46 cm, which corresponds to a capture rate C 10 19 s −1 [7,98]. This is sufficiently fast to ensure that equilibrium will be reached, which implies a DM annihilation rate Γ ann = 1 2 C . If twin pions produced in the dark shower have a typical lifetime τ 1 s, most of them escape from the solar interior before decaying to ∼ GeV diphotons that may be seen by terrestrial experiments. For m π 200 MeV, the incident photon flux is approximately of order 10 −8 cm −2 s −1 , which is comparable to the γ-ray flux seen by the Fermi -LAT quiescent sun data [99], or the flux bound recast from a leptonic final state analysis in Ref. [100]. It would therefore be interesting to examine the future sensitivity of solar observations to this DM annihilation signal. We leave a careful study of these solar constraints to future work.

B. Galactic center constraints
Treating the observed galactic center γ-ray spectrum as an upper bound on the total photon flux due to DM annihilation plus backgrounds, we may instead construct constraints on the m DM -Λ qcd parameter space along the DM relic density theory contour, Ω ψ h 2 = 0.12.
Here we consider a DM annihilation photon spectrum to be consistent with the galactic center data, if the spectrum does not exceed in any bin the total observed photon flux [38] anywhere in its energy range. Note again that m π or f /v is fully determined along the Ωh 2 contour for each point in m DM -Λ qcd parameter space via eqs. (3.3).
The galactic center constraints are shown in Fig. 11. Regions of parameter space previously identified in Fig. 8 as consistent with the GCE spectrum, including the best fit region, obviously still remain viable. Regions with twin pion mass and confinement scales larger than in the best fit region of Fig. 8 feature twin hadron masses that are closer to the DM masses, and therefore their dark showers exhibit a lower multiplicity of hadrons, and hence smaller photon fluxes. In both the Dirac and Majorana cases, such regions remain consistent with the observed γ-ray spectra. For the purpose of setting constraints, the regions for which modeling of the twin showering spectrum becomes unreliable -the red shaded regions -are less consequential than they were in the prior case of carrying out a GCE fit. This is because in regions of parameter space where the γ-ray spectra are significantly below current sensitivities, even substantial variations in the showering spectrum cannot lead to detectable effects. Further, unless the variations in such spectra push them to the point of detectability, these variations can safely be ignored. As a result, the red shaded regions are smaller here compared to those in Figs 8 and 10.
Besides bounds from galactic center data, there are also constraints from the photon fluxes generated by dwarf spheroidal galaxies (dSphs). Some studies (see e.g. Ref. [49]) have found significant constraints on various 2-2 DM to SM annihilation channels, that may otherwise explain the GCE. These results, however, depend on estimates for the astrophysical J-factors of the relevant DSphs, which require confirmation from observation. Future measurements of the kinematics of the member stars in dSphs can resolve this issue, and may provide additional constraints on the dark showering scenario, beyond those shown in Fig. 11.
Finally, let us consider the constraints from γ-ray spectra in the case that m π > 3m π and less than 2m K . In this regime, as discussed in Sec. IV B above, the dominant twin pion decays are to π + π − π 0 or 3π 0 final states, such that a higher multiplicity of softer photons are produced by the dark showers, along with soft muons and neutrinos: The photon to lepton production ratio from twin pion decays is expected to be approximately 11 : 4 in favor of photons. Since six SM final states are produced per twin pion decay in this regime, the mean photon energy will scale down by a factor of three, but the multiplicity will increase by only slightly more than a factor of two. Moreover, any photon produced in a multi-body twin pion decay will have an energy of at most m π 0 /2 in the twin pion rest frame, and therefore the photons produced by such decays cannot be harder, decay-by-decay, than for π 0 → 2γ decays. The astrophysical γ-ray constraints typically loosen as a superlinear power law as the photon energy decreases, but only tighten linearly with multiplicity. Hence in the 3m π < m π < 2m K regime, for which the constraints in Sec. IV D do not apply, one expects the galactic center constraints in Fig. 11 to be universally weaker over the entire parameter space, yielding larger allowed regions.

VI. CONCLUSIONS
We have studied a twin WIMP (T-WIMP) scenario whose twin sector contains a mirror copy of the SM hadrons -i.e. all three twin quark generations with three light quark flavors -but no dark radiation -i.e. no light twin leptons or photons. By comparison with existing twin Higgs frameworks, this scenario represents a relatively unexplored region of the twin Higgs theory space.
Unlike various other twin Higgs frameworks, here the lightest twin degrees of freedom of this hadrosymmetric twin Higgs model are twin pions -pseudoscalars -whose coupling to the SM sector via the Higgs portal is heavily suppressed. While the twin hadron sector can nevertheless be accessed at the LHC through the Higgs portal, any energy injected into the twin sector will just produce a shower of detector-stable, invisible particles. The only robust handles at the LHC, then, are evolving constraints on the invisible width of the Higgs.
From a cosmological and astrophysical point of view, however, the presence of light twin pions results in a much richer set of constraints and detection modes. To avoid a matterdominated or overclosed universe post-BBN arising from metastable twin pions, the decay of the lightest hadron -the π 0 -must be completed by the BBN epoch. For m π 0 < 3m π 0 , this implies the presence of a new SM-twin mediation portal with a mass scale around a TeV or below. This new mediator can be probed at the intensity frontier or at the LHC in almost all of its viable parameter space. Despite the absence of dark radiation and LHC-accessible twin hadron decays, we nevertheless find a remarkable experimental complementarity between the cosmological bounds, the LHC and the intensity frontier.
In the Dirac and Majorana DM scenarios we consider, the dark matter freezes out through twin-electroweak mediated annihilation into twin quarks. Annihilations of T-WIMPs then produce dark showers of twin hadrons which are either stable -such as the π ± and twin proton -or decay to SM photons -the π 0 . Apart from Higgs invisible width bounds, these DM annihilations at the decoupling epoch are competitively constrained by CMB reionization bounds. Such DM annihilations in the galactic center, however, naturally evade all astrophysical cosmic-ray constraints on antiproton or position production and proceed at a rate commensurate with a thermal relic cross-section, that in turn corresponds to a photon flux near the current sensitivities of Fermi -LAT. Simulations of this dark showering process reveal that for the Dirac DM scenario, regions of parameter space that produce the observed DM abundances, as determined within the twin Higgs framework, and as preferred by naturalness and Higgs coupling constraints, are precisely those regions that successfully reproduce the claimed galactic center γ-ray excess in the Fermi -LAT data. Conversely, current galactic center data, when applied as an upper bound on the total γ-ray flux from the galactic center, only partially constrains the available parameter space.
In this work we have mostly considered scenarios that do not admit a significant leptonic decay mode, π 0 → µ + µ − . Should this decay channel be available, the required SM-twin portal can easily be out of reach of any near future experiment. However, the astrophysical bounds or reach for DM annihilations in the galactic center into high multiplicity soft leptonic final states are currently unknown. This may be explored in the future by conducting dedicated dark showering and galactic propagation simulations for this soft leptonic scenario.