Complementarity for Dark Sector Bound States

We explore the possibility that bound states involving dark matter particles could be detected by resonance searches at the LHC, and the generic implications of such scenarios for indirect and direct detection. We demonstrate that resonance searches are complementary to mono-jet searches and can probe dark matter masses above 1 TeV with current LHC data. We argue that this parameter regime, where the bound-state resonance channel is the most sensitive probe of the dark sector, arises most naturally in the context of non-trivial dark sectors with large couplings, nearly-degenerate dark-matter-like states, and multiple force carriers. The presence of bound states detectable by the LHC implies a minimal Sommerfeld enhancement that is appreciable, and potentially also radiative bound state formation in the Galactic halo, leading to large signals in indirect searches. We calculate these complementary constraints, which favor either models where the bound-state-forming dark matter constitutes a small fraction of the total density, or models where the late-time annihilation is suppressed at low velocities or late times. We present concrete examples of models that satisfy all these constraints and where the LHC resonance search is the most sensitive probe of the dark sector.


I. INTRODUCTION
The existence of dark matter (DM) is well-established by observations of its gravitational effects. However, the particle nature of DM is still very much a mystery, despite the ongoing efforts of many complementary experimental searches. Constraints set by XENON [1], LUX [2] and PandaX [3] have strongly ruled out generic DM candidates that interact in a spin-independent manner through a Z-exchange, and are now starting to probe Higgs-mediated interactions (e.g. [4]). These direct detection experiments are complemented by dark sector searches at colliders. The main DM search strategy at the Large Hadron Collider (LHC) is based on missing transverse momentum (MET) balanced by a jet, electroweak (EW) gauge boson or Higgs, known generically as mono-X searches. Searches for dijet or dilepton resonances, while not directly probing the existence of DM, can also effectively constrain models where a mediator particle is responsible for interactions between the Standard Model (SM) and a "dark sector" containing the DM, limiting the parameter space for the mediator. Finally, indirect searches for DM annihilation or decay to SM particles, as well as the well-measured relic abundance of the DM, set powerful limits on the strength and nature of the interaction of DM with the SM. Any model of DM must successfully contend with all of these constraints.
With no hint yet of what the dark sector may look like, we might look to the SM for clues as to its possible composition and structure. In this light, we should * Electronic address: gelor@uw.edu † Electronic address: hongwan@mit.edu ‡ Electronic address: tslatyer@mit.edu § Electronic address: soreqy@mit.edu not be surprised to find bound states in the dark sector; after all, bound states are ubiquitous in the SM, and even the simplest dark sector models with a DM candidate and a force carrier can potentially support the existence of bound states. Dark sector bound states, much like QCD bound states, may be produced when a pair of heavy dark sector particles are produced close to their kinematic threshold and have a sufficiently strong attractive interaction between them. The subsequent decay of these bound states into lighter SM particles can lead to distinctive signatures at the LHC. This strategy has been studied in the context of bound states formed by supersymmetric (SUSY) particles, and has been shown to be a potential search channel at the LHC [5][6][7][8], capable of probing regions of parameter space where traditional searches are challenging.
Dark sector bound states and their potential collider signatures have been studied extensively in the literature. Bound states formed from weakly-interacting massive particles (WIMPs) that are charged under the SM SU(2) L × U(1) Y gauge group or non-SM forces, known as WIMPonium [9], can be detected at the LHC through resonant decays into a pair of leptons, provided the coupling to the mediator which supports the bound state is large enough. Other model-specific dark sector bound state collider searches that have been proposed include searches for higgsino bound states in λ-SUSY and bound states within the self-interacting DM framework [10]; DM bound states in a U(1) vector portal model decaying into multilepton final states, which can be searched for at B-factories [11]; and a Higgs portal model with decays to electrons which can be searched for at the LHC [12]. Mono-photon searches at lepton colliders can also potentially be used to probe the full resonance structure of the dark sector [13]. However, the large couplings typically required for detectable bound states often predict large signals in direct detection experiments, espe-cially if the light force carrier responsible for the bound state formation also couples to the SM; likewise, in this light-mediator regime, searches for the mediator are often a more promising dark-sector discovery channel than searches for the bound states [14].
In this paper, we broadly explore the challenges of building a dark sector model which can be discovered through the production of a bound state at the LHC, in light of the current stringent and complementary experimental constraints. Direct detection limits can be evaded in models with TeV-mass DM if the DM candidate only has an off-diagonal coupling to the SM that couples the DM, the mediator and a heavier dark sector state, so that at tree-level, the DM only scatters into this heavier state when interacting with the SM [15,16]. At the LHC, dark sector particles can be produced on their kinematic threshold and form a bound state B, which can subsequently undergo annihilation decay into a pair of SM leptons, showing up as a dilepton resonance at the LHC. 1 We will show that in models where the mediator between the SM and the dark sector couples to two different states in the dark sector, it is possible to arrange for such a resonance to occur and have a substantial branching ratio into SM leptons. In these scenarios, searches for a dilepton resonance from B are complementary to the existing mono-X and vector resonance searches that are already deployed for dark sector searches at the LHC, with the ability to probe higher mass scales for the mediator and DM. Since B can have the same quantum numbers as the SM mediator, we explore the importance of mixing between bound states and mediator particles with equal quantum numbers and similar masses.
Models with bound states that are detectable at the LHC can also possess large indirect signals, as the longrange potential implied by the existence of bound states generically enhances the annihilation cross section for slow-moving DM particles, and the bound state formation and decay can also serve as an annihilation channel. We will study the constraints from indirect detection and cosmology that result from considering these effects.
The rest of the paper is structured as follows. In Sec. II, we will make some remarks on the general features of dark sector models where the bound-state resonance search is the most sensitive channel. We will discuss why the bound-state resonance search is complementary to the current dark sector search strategies used by the LHC experiments for such models, and discuss their general phenomenology in direct, indirect and collider DM searches. In Sec. III we will lay out some specific models containing bound states in the dark sector and study their phenomenology. We will first discuss the MSSM in the pure wino/higgsino limit, which already meets some 1 Di-jets are also a plausible search strategy, but the backgrounds and triggers make this much more challenging to explore.
of the criteria needed for a successful model with bound states, although the production rate at the 13 TeV LHC is too small for detection. We will then discuss two vector portal models which realize the requirements needed for a viable dark sector with bound states to be probed by the LHC. In Sec. IV we compute and discuss the potential experimental signatures of these models. Our conclusion will then follow in Sec. V.

II. PHENOMENOLOGY OF BOUND STATES
The existence of DM bound states has implications for the phenomenology of the dark sector, and for its signatures in direct, indirect and collider searches. In this section, we consider the circumstances under which collider searches for bound states can probe otherwise unexplored regions of DM parameter space. Aside from these searches, DM bound states with long lifetimes have also recently been shown to have potentially interesting implications for neutrino experiments [17].
As we will show, models where bound-state resonance searches at the LHC probe new regions of parameter space are most easily realized in the presence of several common features: 1. DM couples to at least two distinct force carriers; one of these, Y , is light and mediates the bound state formation, while the other, V , is heavier and couples appreciably to the SM. The constraints from LHC resonance searches of the bound state are most competitive when the SM mediator V is heavier than twice the DM mass; subsequently undergoing a similar decay. The coupling of the mediator between the dark sector and the SM to quarks (gq) and to the DM (gχ), as well as the coupling responsible for the Yukawa potential that forms the bound state B (αB ≡ g 2 B /4π) are shown. In our models, V is always a vector, while Y can be either a scalar or a vector.
we will denote the heavy mediator as V and its mass by m V , and the light mediator by Y and its mass as m Y (for "Yukawa"). In the example models we present, V will be a vector in all cases, but Y can be either a scalar or vector. In principle, V could also be a scalar (or a scalar bound state can mix directly with the Higgs sector [12]), but we will leave the analysis of such scenarios to future work; as we will see, a vector mediator facilitates a sizable production cross section and a large branching ratio to leptons, while evading direct detection bounds.
Many of the earlier works in the literature on bound states exhibit some of these features. Both [9] and [10] introduce an additional mediator to support the bound state formed from DM charged under the EW gauge group, so that the couplings between the DM to the light mediator can be made large. An additional mediator was also introduced in [17] to alleviate the tension between a suitably light mediator that can support a bound state and the need for a massive enough SM mediator that can decay into electron pairs. In [11], direct detection limits are avoided by having sub-GeV DM. Furthermore, there is only one vector boson to mediate both the bound state formation and the interaction with the SM, at the cost of allowing the bound state to decay into 4-or 6-lepton final states. This is an important signature in B-factories for DM with a mass on the order of a GeV [11]. In principle, this scenario can be probed at the LHC by multi-lepton searches, or by di-photon searches where two e + e − pairs are detected as fake photons [12]. However, multi-lepton signatures turn out to be relatively unimportant for the kinetic mixing models that we will study later.
We will demonstrate that the characteristics listed above can be achieved in Higgsed dark-sector models, with a vector portal between the dark sector and the SM. But before we give examples of such models, we will first discuss each of these criteria in more detail.

A. General Model Building Considerations
The existence of DM bound states in a Yukawa potential with range 1/m Y is only possible if [11,18]: where m χ is the DM mass. Thus, the presence of a bound state supported by scalar or vector exchange requires a relatively light force carrier -certainly lighter than the dark matter itself, for weak couplings. For more complex dark sectors with potentials that couple multiple two-particle states (e.g. the neutralino sector of supersymmetric models), the details of this criterion may be modified, but it is still generically true that there must be a force with range longer than the Bohr radius of the bound state, i.e. there should be at least one mediator with m Y α B m χ . If this force carrier is also the mediator between the DM and the SM, then searches for the force carrier will generally offer a more accessible probe of dark-sector physics than searches for the heavier DM, both because the force carrier is lighter and because it couples directly to SM particles (see e.g. [14]). This leads us to consider models where there are at least two distinct particles that couple to the DM, one which has appreciable interactions with the SM (and can be heavier than the DM itself), and the other of which mediates the bound state formation and so must be light.
One alternative to this structure is the case where the DM is charged under the SM SU(2) L EW gauge group, and the photon, W and/or Z support the bound state; this is possible, for example, for bound states consisting of neutralinos and/or charginos [19]. However, as we will show later in this work, at present-day colliders the production rate for such EW bound states is undetectably low.
Returning to dark-sector models with at least two mediators, the presence of the light mediator has some immediate implications. First, the DM will generically annihilate into the light mediators. If these mediators are absolutely stable, they will constitute some fraction of the DM relic density, which must be sufficiently small; if they are below the ∼ MeV scale in mass, they may be constrained by limits on the number of effective relativistic degrees of freedom in the early universe (e.g. [20,21]). We will generally assume that these mediators decay through some small mixing with the SM, on timescales less than one second, so that they do not affect Big Bang nucleosynthesis; in this case, while the coupling can be made small enough that these mediators do not contribute substantially to collider signals and direct detection, indirect detection constraints from this annihilation channel must be considered.
Resonance searches for bound states will typically become difficult when there is a significant branching ratio of the bound state into dark sector states (while the dark sector states could decay promptly to SM particles as in Ref. [11], for the heavy DM models considered in this paper, which are most relevant for LHC searches, we expect this signature to be relatively unimportant, as we will explain later). Thus, a model where collider searches for bound state resonances are effective must ensure that the branching ratio of the bound state (formed at a collider) into two light mediators is small relative to the decay of the bound state via the off-shell heavy mediator into a pair of SM particles. One way that this can be achieved is if both the heavy and light mediators are vectors, then suppression of the spin-1 bound state decay to two light mediators is automatic: charge parity symmetry forbids the decay of a spin-1 s-wave bound state into two vectors, so any decays into dark sector vectors must be a 3-body process. In fact, decays into any number of the light mediators can be completely forbidden if the bound state is formed not from a particle-antiparticle pair, but from two different fermions in the bound state with nontrivial quantum numbers, which cannot be conserved if the bound state could decay into states containing only light mediators. This behavior is natural in cases where the mediator V couples off-diagonally to the multiplet containing the DM. Models of this type have additional advantages in evading constraints from direct detection.
Note that if the mediator to the SM is a vector, then the bound states formed at the LHC by resonant production will dominantly be spin-1 s-wave states; if the mediator is a scalar, they will instead dominantly be spin-0 s-wave states. The spin and angular momentum of the bound states determine their possible decays.

B. Vector-Bound State Mixing
When V and B have similar masses or the coupling between V and the constituents of B is large, significant mixing can occur between the two states if they have the same quantum numbers. Both the V -resonance and B-resonance diagrams in Fig. 1, together with higher order diagrams with more inter-conversions between V and B, need to be re-summed. The new mass eigenstates that result from the mixing have masses and widths that are shifted with respect to their unmixed values by an amount determined by the strength of the mixing.
The formalism that accounts for the mixing was used to study Z-toponium mixing [22][23][24][25][26], and more recently to study Higgs-stoponium mixing [27]. The mixing shifts the masses and widths of the unmixed states, denoted V 0 and B 0 , to new values given by the eigenvalues of the following mass matrix: where all masses and widths are for the unmixed states, and f is a model-dependent parameter determined by the coupling between V 0 and B 0 .
If f is small compared to the difference in the diagonal entries (see Eq. (4) below), the final mixed states V and B are approximately their respective initial unmixed states, up to higher order corrections. The width of V 0 should be evaluated at the appropriate energy scale √ s at which the final mixed resonances V or B are produced; this scale dependence is important especially when m B lies below the χχ open production threshold while m V,0 lies above it. The width of B 0 should not include decays through mixing with the V : such effects are exactly what the mixing accounts for. For the kinetic mixing models that we will consider later, we take Γ B,0 = 0, since the dark sector particles do not have any tree-level coupling to the SM, and the unmixed width of the bound state excluding mixing into the SM is always much smaller than Γ V,0 . After mixing, the mixed mass eigenstates are rotated by a complex mixing angle θ with respect to the unmixed states, and the masses and widths are shifted by [23,25] Q V = Q V,0 cos 2 θ + Q B,0 sin 2 θ + f sin 2θ, where The rotated mass eigenstate B therefore develops a coupling to the SM through its V 0 component. When the mixed masses m V and m B are nearly equal, a resonance search for each individual mass eigenstate becomes impossible, since the s-channel diagrams with intermediate V -and B-states interfere with each other, and the end result is a cross section that may not have a Breit-Wigner form. However, if θ is small, Eq. (3) shows that the mixed mass eigenstates are separated by ∆m 2 ∼ 4f Re(θ), where Re(θ) is the real part of θ. Furthermore, the shift in the masses defined by Eq. (3) and (4) always results in a mass eigenstate that is lighter than both m V,0 and m B,0 , and is therefore always strictly below the threshold for open production of χχ. These two facts can ensure that the lighter resonance is always narrow, as it cannot decay into χχ, and is always wellseparated from the heavier resonance. We have checked that this is always the case for the models that we consider later.
Finally, in the limit of small θ, this mixing procedure gives a final decay width Γ B that agrees with the perturbative calculation to O(θ 2 ), i.e. with the result obtained by summing the partial widths of B 0 decaying through mixing with V 0 (with Γ V,0 evaluated at s = m 2 B,0 ), which then decays into SM final states [22,23]. Throughout this paper, we will therefore qualitatively discuss the nature of the B resonance using the perturbative picture, while taking the mixing fully into account quantitatively. We will also not make a distinction between B 0 and B or V 0 and V , unless we are explicitly discussing the mixing.

C. Collider Signatures
There are three important classes of collider signatures for models of the type we have discussed: (i) mono-X, where the DM state χ is produced and observed as MET recoiling against a SM final state such as X = j, h, W, Z ; (ii) V resonant production with decaying channels such as dilepton, dijet or any other SM final states, and (iii) B resonant production with m B ≈ 2m χ , decaying into a pair of leptons or jets.
The three channels probe different physics, as well as different regions of the dark sector parameter space. The mono-X channel is an unavoidable signature of DM. The properties of the B resonance are completely determined by the DM mass and its self-interaction through the light mediator Y ; therefore, by analyzing its properties, we study the DM directly. The V resonance on the other hand probes the structure of the dark sector, but is not directly related to the puzzle of DM.
The mono-X signature has been discussed previously [28][29][30][31][32], and there are on-going searches at the LHC. We will demonstrate that for the models we consider, mono-jet searches probe a different region of parameter space than bound state resonance searches.
The production rate for the bound state B at a pp collider is given by (see [7,10,33] and also Appendix A) where J is the spin of the bound state; for a bound state produced from a vector mediator, J = 1. ζ(s) is the Riemann ζ-function, which takes into account the cross section for the production of all of the excited states of the bound state. Here τ B ≡ m 2 B /s, m B is the mass of the bound state, and √ s is the collider center-of-mass energy.
L qq is the parton luminosity function defined as with f q (x) being the parton distribution functions (PDF), taken from [34] for calculations in this paper.
In the perturbative limit, we can write where g q (g χ ) sets the coupling of the mediator V to quarks (DM) and ψ(0) is the wave function of the bound state at the origin. For a Coulomb-like potential with coupling α B (i.e. where the mass of the bound state mediator m Y can be neglected), |ψ(0)| 2 = α 3 B m 3 χ /8π. This perturbative B production cross section can be understood in three limits: (i) the heavy mediator limit, m V m B ; (ii) the light mediator limit, m V m B , and (iii) m V ≈ m B . The cross section in each limit is These equations show that the B production cross section is enhanced when its mass is close to the mediator mass, and suppressed in the other two limits. Thus, we expect stronger sensitivity in this channel when m V ≈ m B . Moreover, if m B m V , which is in the limit where V can also support dark matter bound states, the B production cross section is suppressed by Γ 2 V /m 2 B relative to the m V ≈ m B region. We also can see that for models where B is heavy enough to decay primarily into two or three V 's which then decay into 4 or 6 leptons at the LHC, the production rate of the bound state is suppressed relative to the regime where m V ≈ m B . The mediator production cross section, V , is where τ V ≡ m 2 V /s. As we pointed out above, the V resonance search does not directly probe the dark matter content. Further searches must be used to uncover the dark sector after discovering the mediator between the SM and the dark sector. Most importantly, when m V > 2m χ and g χ g q , g (the coupling to leptons), the branching ratio of V to SM particles becomes small, and resonance searches for V grow ineffective. The full mixing calculation also bears out this conclusion: once m V,0 > 2m χ , the V resonance is heavier and lies above the χχ threshold and is a wide resonance, while the lighter B resonance remains narrow and below the threshold.
The comparison between mono-X and bound state production is more complicated as the backgrounds for the two searches are different, and a more detailed comparison is required; we will show results for some specific models below. On generic grounds, the mono-X cross section is reduced because of the PDF price of the additional jet. However, the two production cross sections scale as α s g 2 q g 2 χ and α 3 B g 2 q g 2 χ , for the mono-X and bound state cases respectively. Thus for α 3 B α s we expect a reduced sensitivity in the bound-state searches; this suggests α B rather close to 1 will be required to make bound-state searches competitive. Moreover, the monojet search becomes ineffective once m V < 2m χ , since the mono-jet process must then proceed through an off-shell V .
In summary, the mono-jet search probes the region of parameter space where m V > 2m χ , while the V resonance search is more sensitive to the region where m V < 2m χ . The bound-state production cross section, on the other hand, is enhanced precisely in the intermediate region, and outperforms the other two searches when m V 2m χ . These three searches are thus complementary, and probe different parts of parameter space, as we will show explicitly in our models below.

D. Direct Detection Limits
Direct detection searches are very sensitive probes of DM, especially for DM with substantial couplings to hadrons, and mass at the EW scale or higher. Thus, viable models of dark resonance signals at the LHC must evade direct detection bounds.
A naive estimate of the DM-nucleon scattering cross section at tree level, in terms of the parameters discussed in the previous subsection, gives σ ∼ g 2 q g 2 χ m 2 N /m 4 V ∼ 10 −40 cm 2 g 2 q g 2 χ (TeV/m V ) 4 , assuming m V is much larger than the typical momentum transfer in the scattering, and m χ is much larger than the nucleon mass m N . For comparison, under standard assumptions, the limit from XENON 1T on this scattering cross section is of order 10 −45 cm 2 (m χ /TeV) [1]. Thus, if the elastic scattering spin-independent cross section is unsuppressed, we infer that the product of couplings g 2 q g 2 χ 10 −5 m 4 V m χ /TeV 5 . This simple estimate is broadly consistent with more carefully obtained limits on a dark sector interacting with nucleons through a vector mediator for current and future direct detection experiments [35,36]. Reasonably large couplings and sufficiently low dark sector masses are necessary for the significant production of the bound state resonance, but this parameter region of interest (g q g χ ∼ 1 and m V ∼ 2m χ ∼ 1 − 4 TeV) is generically in tension with direct detection bounds.
However, any suppression to the naive tree-level cross section can alleviate this tension. As mentioned above, a simple scenario ("inelastic dark matter") that leads to suppressed direct detection signals posits that the coupling between the DM χ 1 and the mediator V involves an unstable partner particle χ 2 , and the mass splitting between the DM and its partner is greater than the maximum kinetic energy of DM particles in the halo [15,16]. detection of the DM particles. Such models also have interesting consequences for bound state formation at the LHC: if the bound state is produced in the s-channel from the mediator V , it will automatically be composed of the DM and its partner particle, or may only involve dark sector particles in the same multiplet as the DM, and not the DM at all. In such models, elastic scattering can still occur, but only at loop level. The direct-detection spin-independent cross section for scattering off a nucleon with target mass m T is given by [37] where n p,n are the number of protons and neutrons respectively, and f p,n are the corresponding matrix elements. We can generalize the effective operator analysis of [38], to the dark sector models that will be of interest to us. Then where N = n, p, and we sum over u, d, s quarks. Here the first term comes from a 1-loop diagram involving the Higgs, while the second term is a box diagram with two V propagators. q(2) andq(2) are the second moments of the quark and anti-quarks PDFs and f Tq = N |m qq q|N /m N is the nuclear form factor. For these we use the numeric values from [38]. In the wino scenario discussed in [38] these two contributions f q and g 1,2 q are non-negligible but of opposite sign, thus leading to a cancellation. In the case of the dark sector models that we will consider later, the second contribution is suppressed by the small coupling between the SM mediator and the SM, while the first contribution will be negligible as the dark sector coupling to the SM Higgs will always be very small. Explicitly, we can write g (1,2) are functions computed in [38], is the small mixing parameter, and Q is the charge of the quark. We find that in regions of parameter space of interest to the present work the contribution from loops to direct detection is thus no larger than σ SI ∼ 10 −48 cm 2 and is thus unconstrained.

E. Overclosure and Indirect Searches
In general, if annihilation to the light mediators that support the bound state is allowed, this process will tend to dominate freezeout. The same attractive potential that permits bound state formation will also generically enhance annihilation through the Sommerfeld enhancement [39,40], potentially giving rise to large indirect signals in the present day. Formation of bound states followed by their decay can also significantly enhance indirect signals (e.g. [41]), if the mediator supporting the bound state is light enough that radiative capture of two DM particles into the bound state is kinematically allowed.
Let us first note that there are several possible annihilation channels which are p-wave suppressed at late times [42]; if these processes dominate freezeout, the latetime indirect detection signals will generally be very suppressed. We will see an example of this when we consider a model where the dominant annihilation is of Majorana fermions to light scalars. Furthermore, if DM-DM scatterings experience a repulsive potential rather than an attractive one, the DM-DM annihilation will be exponentially suppressed at low velocities [43]. Note that in order for bound-state searches to be interesting in such a scenario, there must be at least one other particle with which the DM can form a bound state, and the DM must have an attractive interaction with that particle. In this initial discussion, therefore, we assume the dominant annihilation is s-wave and experiences an attractive Sommerfeld enhancement, to explore when indirect searches can set interesting constraints.
Let us consider the simple case where the effective potential experienced by the DM is a Yukawa potential, as discussed above. The s-wave Sommerfeld enhancement for a Yukawa potential with coupling α B and mediator mass m Y can be well approximated by [44,45]: where v rel is the relative velocity between the DM particles, and in the second line we use cos θ > −1 with an angle controlling the resonance positions. The inequality is saturated for real values of θ, at the minima between resonances where cos θ = −1. (It is also approximately saturated where m Y → 0 and v rel α B , where S ≈ 2πα B /v rel .) Note that for fixed m χ , requiring the correct relic density fixes α B , if the assumptions are made that (a) this channel dominates during freezeout, and (b) the mass of the mediator is irrelevant during freezeout. The latter assumption is approximately true away from resonances and if m V /m χ is smaller than the typical velocity of particles around freezeout (v ∼ 1/3). For large values of m Y (requiring large α B , since m Y /m χ < α B /1.68), or values of m Y corresponding to resonant Sommerfeld enhancement (cos θ → 1 as v rel → 0), freezeout may be more complicated and needs to be studied more carefully; we will include the full m Y dependence when we examine specific models.
However, if we consider α B to be fixed given m χ , and hold v rel fixed, then our expression for the lower bound on the Sommerfeld enhancement is a monotonically decreasing function of m Y ; thus indirect detection will set a lower bound on m Y (all values of m Y below this threshold will be ruled out). Since the requirement for bound state formation sets an upper bound on m Y , one can ask whether these two criteria are in conflict.
In order for the model to avoid exclusion by indirect detection (except possibly where m Y is important to freezeout), this minimal Sommerfeld enhancement must be permitted by the data. Note that for α B v rel ( v rel ∼ 10 −3 in the present-day Milky Way halo), this minimal enhancement will reduce to an α B -and v rel -independent prefactor of: This minimal Sommerfeld-enhanced cross section is rather close to indirect detection bounds for a DM species that comprises 100% of the DM and whose abundance is set by thermal freezeout, for DM masses below ∼ 1 TeV (e.g. [46,47]); permitting m Y m χ v rel would generally significantly overproduce limits from indirect detection, unless Y decays primarily into invisible channels. If we assume m Y m χ v rel in the present day, then we can approximate S 6α B m χ /m Y , and thus if the maximum allowed Sommerfeld factor is S max , then m Y 6α B m χ /S max . Of course, smaller values for m Y are permissible if the species that forms bound states comprises only a small fraction of the overall dark matter density.
If the dominant annihilation channel consists of swave annihilation to mediators coupled to the DM with strength α B , then the annihilation cross section at low velocities is of order σv rel ≈ πα 2 B /m 2 χ (this expression is exact for Dirac or pseudo-Dirac DM annihilating to U(1) dark gauge bosons). Requiring that this cross section fall below the thermal value of σv rel ≈ 2 × 10 −26 cm 3 /s ≈ 1.7×10 −9 GeV −2 suggests an overclosure bound of m χ α B ×43 TeV. As we will see, we will generally be interested in masses around a few TeV and α B 0.1, so the overclosure bound will not typically be particularly constraining. This estimate ignores Sommerfeld enhancement and bound state formation during freezeout, which can be important [48,49]. For α B 0.1, the Sommerfeld enhancement is non-negligible during the freezeout epoch; however, for attractive Sommerfeld enhancement, including this effect only reduces the late-time relic abundance. This further relaxes the overclosure bound, and since it reduces the abundance of the species in question, also weakens constraints from indirect detection. (However, it makes it more challenging to generate 100% of the DM abundance by the same species that forms bound states.) Likewise, radiative formation of bound states can also contribute to the depletion of DM at early times and indirect signals at late times [41,[49][50][51]. These radiative processes are only kinematically unsuppressed if enough energy is available to produce an on-shell light mediator, i.e. the binding energy + kinetic energy of the particles is greater than m Y . Bound state formation can also occur through radiation of an off-shell heavy mediator that decays to SM particles, but such processes will be suppressed by a small mixing with the SM and also by the mass of the heavy mediator. Thus, there are two distinct regimes for m Y from an indirect-detection perspective: , where bound states exist but radiative capture into them is suppressed, and m Y α 2 B m χ /4, where radiative capture processes are unsuppressed. We will ignore bound-state effects in the former case, but account for their impact on indirectdetection signatures in the latter case.
However, we will ignore the effects of bound-state formation during freezeout. A careful treatment of boundstate effects during freezeout requires accounting for dissociation of the bound states through interactions with the light-mediator bath. If m Y α 2 B m χ /4, then for α B 0.5 we expect the temperature at freezeout to be comparable to or larger than the binding energy (taking the standard estimate T freezeout ∼ m χ /20), and so dissociation effects could be substantial. Thus while the presence of radiative capture into bound states during freezeout may further deplete the DM abundance, relaxing both the overclosure and indirect limits further, a full calculation would require a careful analysis (as performed in e.g. Ref. [49]).
We will show that the indirect detection constraints and overclosure limit cannot fully exclude the regions of parameter space relevant to collider searches for the bound states, even without taking the impact of boundstate effects on freezeout into account, for both models we consider. Since including the bound-state effects during freezeout would only relax these constraints further, we are justified in neglecting them for purposes of this work.

F. Dark Matter Self-Interactions
Constraints on DM self-interactions require that σ/m χ 1 cm 2 /g ≈ 1/(60MeV) 3 [52][53][54]. As we are interested in the regime where a long-range potential exists and can support bound states, we cannot use the Born approximation to estimate scattering rates. However, if m Y /m χ v rel while still satisfying Eq. (1), the typical relative velocity of DM particles in galaxies and galaxy clusters, then we can make the approximation that swave scattering dominates and use the analytic estimate for the scattering cross section derived in Refs. [55,56].
The scattering cross section is approximated in the low-velocity limit by [55,56] where We see that away from resonances, which occur when cot(π √ c) diverges, the size of the phase shift . The regime where the s-wave contribution dominates is thus a regime where (away from resonances) this phase shift is small, and we can write: which is just the geometric cross section. Assuming this geometric cross section, we see that the self-interaction bound will be satisfied provided (m χ m 2 Y ) 1/3 100 MeV, which for 1 TeV DM requires only that m Y 1 MeV. Thus, away from points in the parameter space where there is a near-zero-energy bound state, we expect the self-interaction rate to be undetectable, despite the rather large couplings we invoke.

III. DARK SECTOR MODELS WITH DETECTABLE BOUND STATES
We now consider two examples of phenomenologically viable DM models containing bound states, which can lead to interesting signatures at the LHC. These models serve as examples of how to build non-supersymmetric models that realize the requirements laid out in Sec. II. We will show that in these models, the search for boundstate resonances can probe parameter space which is not accessible to mono-jet searches and resonance searches for the mediator.
In the first model, which we label as the "pseudo-Dirac" model, the dark sector consists of a pair of almostdegenerate Weyl fermions that are charged under a darksector U(1) D gauge group, which is broken by a dark Higgs-like scalar. These fermions can form bound states with the dark Higgs as the mediator. The second model, which we refer to as the "triple Higgs" model, is based on a completely broken SU(3) D gauge theory, with the dark matter candidate being a Dirac fermion in the fundamental of the gauge group. Much of the phenomenology of this model, including bound state formation and couplings to the SM, is derived from the symmetry breaking pattern of the theory, with both the mediator that supports the bound state and the mediator to the SM being massive gauge bosons of the broken SU(3) D group. In both cases, the dark sector interacts with the SM via a vector portal with kinetic mixing, and the DM direct detection cross section is suppressed by the fact that at tree-level the DM scatters into a heavier state.
Before introducing these models, however, we will consider a simpler scenario that is familiar from SUSY, that of pure wino/higgsino DM (Sec. III A). We will show that the production rate of wino/higgsino-onium bound states at the LHC is too small to be constraining, but this scenario shares many of the properties of our more complicated dark-sector models and thus has pedagogical value. We will then review the details of kinetic mixing between new dark gauge bosons and the SM neutral gauge bosons (Sec. III B), since this mechanism describes the leading interaction of the SM with the dark sector in both darksector scenarios we consider, before describing in detail the two models (Sec. III C and III D).

A. A Weakly Interacting Example: SU(2)L Minimal Dark Matter
Sub-TeV superpartners of the EW bosons and of the two Higgs doublets in SUSY theories can potentially be produced and detected at the LHC, with the lightest neutralino being a particularly well-motivated, weakly-interacting DM candidate. Outside of SUSY theories, models of "minimal dark matter" where the DM transforms under a low-dimensional representation of SU(2) L have similar phenomenology to neutralino DM [43,57,58]. Pure wino or higgsino DM corresponds to the lowest-lying mass eigenstates from, respectively, an SU(2) L triplet of Majorana fermions or an SU(2) L Dirac fermion doublet with hypercharge 1/2. The hyperchargezero SU(2) L quintuplet is also a viable "minimal dark matter" candidate.
If the DM transforms as part of a SU(2) L multiplet, then it will be accompanied by heavier charged partner particles in the same multiplet. After EW symmetry breaking, the wino triplet separates into a lighter neutral Majorana fermion χ 0 and a heavier charged Dirac fermion χ ± ; the higgsino multiplet gives rise to two neutral Majorana states χ 1 , χ 2 and a charged Dirac fermion χ ± . These charged partners can always form Coulombic bound states; when the DM is sufficiently heavy, W and Z exchange may also support bound states including the DM itself (e.g. [19]). Numerical calculations indicate that for wino DM there is a crossover point at a DM mass of around 5 TeV, where the ground state transitions from being primarily composed of χ + χ − bound by photon exchange, to being composed of an admixture of χ 0 χ 0 and χ + χ − bound by the gauge bosons of an approximately unbroken SU(2) L symmetry. SU(2) L -DM models have many attractive features of the type discussed in Sec. II, and behave as prototypes for the models of interest to us. They naturally possess multiple mediators, one of which is massless and sup-ports bound states, while the other massive mediators are all known particles in the SM. The SU(2) L multiplets contain several states nearly-degenerate with the DM; the couplings of the gauge bosons with the DM and its partners are naturally off-diagonal, and so the elastic scattering relevant to direct detection proceeds only at one-loop level (and also suffers from additional cancellations which suppress the rate further [59]). Direct and indirect constraints on wino and higgsino DM have been studied extensively; thermal wino DM constituting 100% of the DM is in tension with H.E.S.S. observations of the Galactic Center (e.g. [60][61][62][63][64][65][66]), but a subdominant wino DM contribution at lower wino masses is difficult to exclude. Pure higgsino DM is not currently experimentally testable by either direct or indirect detection [67,68].
A complicating factor in SU(2) L DM models is the presence of multiple mediators that can potentially support a bound state, which become most important if the DM is heavy enough that α W m χ m W , m Z with α W = g 2 W /4π ≈ 1/30. In this case, there is a longrange potential that mixes the two-body DM-DM state with other particle anti-particle states, i.e. χ 0 χ 0 mixes with χ + χ − in the wino case, and χ 1 χ 1 can mix with χ 2 χ 2 and χ + χ − in the higgsino case. This can lead to χ + χ − states that are only pseudo-bound, despite the presence of the photon-mediated Coulomb potential: if the combined W/Z/γ-exchange potential is not deep enough to also bind the χ 0 χ 0 component, then the χ + χ − state (or e.g. χ ++ χ −− in representations, such as the quintuplet, where higher-charge states exist) may decay rapidly to unbound χ 0 χ 0 through t-channel exchange of W bosons. Parametrically, the cross section for where ∆ is the available energy (i.e. the splitting between the mass of the χ + χ − two-body state, including any binding energy, and the mass of the final χ 0 χ 0 state). By comparison, the cross section for annihilation to SM quarks, leptons and gauge bosons is of order σv rel ∼ α 2 W /m 2 χ . Thus we expect the former to dominate over the latter However, there is an important caveat to this argument: in fermionic models of this type, this mixing between χ 0 χ 0 and χ + χ − occurs only in the states with even L + S (where L and S are the quantum numbers describing the total orbital angular momentum and total spin of the bound state); states with odd L + S have symmetric wavefunctions and cannot support two identical fermions. Since the mediator to the DM is an EW gauge boson, the bound state dominantly produced at colliders has L = 0 and S = 1; these are true bound states, not pseudo-bound, and cannot decay rapidly to pairs of identical DM particles. In particular, in the pure wino case the χ + χ − bound state with L = 0, S = 1, denoted B w , decays dominantly via an s-channel γ/Z to SM fermion pairs or through a t-channel exchange of a χ 0 into a W + W − final state [19]; final states involving the DM are suppressed. We will see this behavior arise again in our example dark-sector models.
The pure higgsino limit serves as an example of a model where there are two neutral mass eigenstates that can be close in mass, denoted as χ 0 1 and χ 0 2 , the lighter of which (χ 0 1 ) is the DM. In this case, the decay of χ + χ − to χ 0 1 χ 0 2 may be allowed. If ∆ +0 ≡ 2m χ ± − m χ 0 1 − m χ 0 2 < 0, the χ + χ − bound state never mixes into the χ 0 1 χ 0 2 from kinematic considerations. When ∆ +0 > 0 however, the χ + χ − can simply decay into free χ 0 1 χ 0 2 , and if the width for this decay is significantly larger than the width of the χ + χ − bound state, the bound state is effectively never formed. 2 Thus, for the pure higgsino case the sign of the parameter ∆ +0 is critical to the bound state phenomenology, at least for DM masses below the TeV scale. This parameter is positive when the lightest neutralino is a pure higgsino and both the wino and the bino are taken to be infinitely massive, but there exists a range of SUSY-breaking parameters which can produce a lightest neutralino that is almost purely higgsino with a significantly more massive bino and wino, while having ∆ +0 < 0 [69]. With this choice, a χ + χ − higgsino bound state B h can be formed and can decay in the same way as B w , albeit with different coupling constants to the EW bosons.
Unfortunately, if the DM is part of an SU(2) L doublet or triplet, the bound state production rate at the LHC is too small to be observed. This is due to the smallness of the EW couplings, which controls the production rate. Figure 3 shows the production cross section times branching ratio into leptons of chargino-onium states for fermions charged under the EW gauge group in different representations. Chargino-onia from both pure winos and pure higgsinos have production cross sections that are far too small for dilepton searches at the LHC to be effective. However, for DM in a larger representation of SU(2) L , fermions having large electromagnetic charges Q can be produced. The production cross section of these states scales rapidly with Q, while the partial widths into SM particles remain unchanged. The enhancement factor relative to the pure wino is Q 8 , with Q 6 coming from the wavefunction of the bound state at the origin |ψ(0)| 2 , and an additional Q 2 from the coupling of these fermions to the γ and Z. For charginos with Q = 4 in an SU(2) L 9-plet, the production cross section for the χ 4+ χ 4− chargino-onium becomes large enough to be probed by the current dilepton resonance search results. Such large representations are generally disfavored since they lead to non-perturbative values of α W below the Planck scale [57]; however, these results more broadly demonstrate that models with large coupling constants or large charges are particularly suited for bound state searches at the LHC. Searches for multi-charged lepton bound states decaying into two photons, for example, 2 When the widths are comparable, bound state decays into χ 0 1 χ 0 2 becomes an additional decay channel, together with decays to W + W − or SM fermions. have been shown to be effective in searches for leptons with a sufficiently large hypercharge [70].

B. Kinetic Mixing
We now turn our attention the dark sector models that we briefly described above. Both of the models we will consider interact with the SM through a vector portal, with kinetic mixing with the SM U(1) Y : where V µν (B µν ) is the field strength of the dark gauge boson (SM hypercharge), and we have included the mass term for both V and the SM Z. Here, V µν can be nonabelian: such a mixing term appears in the triple Higgs model in the form of a dimension-5 operator H a D V a µν B µν where H a D is an adjoint scalar that acquires a VEV, and a = 1, · · · , 8 is an SU(3) D color index.
This interaction can be diagonalized in the mass basis; a detailed description of this diagonalization procedure is discussed in [71,72], and reviewed in Appendix B. In the non-abelian case, only the abelian portion of the field strength is diagonalized, with the non-abelian portion remaining as an interaction term in the model. The diagonalization introduces an -suppressed coupling between the physical dark gauge boson and the SM electromagnetic, J µ EM , and weak-neutral, J µ Z , currents, as well as an -suppressed coupling between the SM Z-boson and the dark sector current, J µ D : where A is the SM photon, s W (c W ) is the sine (cosine) of the weak mixing angle, and r ≡ m Z /m V . All of the fields are given in the mass basis: note that the DM fermionic current couples directly to the Z, so both V and Z mediate the production of dark sector particles with qq interactions, and both must be included in amplitude calculations.
The mixing between V and Z also shifts their masses by a fraction of O( 2 ): the shift in the Z-mass has important consequences for EW precision constraints on these models which we will discuss below, but otherwise these shifts will be neglected for the rest of the paper. We will always assume that r 1 throughout in both models.

C. U(1)D Pseudo-Dirac Dark Matter
We now consider a simple, viable dark matter model, where the bound state signature gives complementary information about the dark sector and probes different region of the parameter space than the mono-X searches. Our model is based on the "minimal model" of [73] (loosely based on the "excited dark matter" scenario of [74]), but we use an ordinary Yukawa interaction between the dark Higgs and the fermions in the dark sector instead of a dimension-5 operator.
This model contains a gauged U(1) D field, V , kinetically mixed with the SM U(1) Y , a Dirac fermion Ψ and a dark Higgs, which in unitary gauge can be written as with v D as its VEV. The U(1) D charges for the fermion Ψ and Φ D are 1 and 2 respectively. The Lagrangian is where D µ ≡ ∂ µ − ig D V µ is the covariant derivative for Ψ and D µ ≡ ∂ µ − 2ig D V µ is the covariant derivative for Φ D , with C denoting charge conjugation. Following [75], we write Ψ as a Weyl fermion pair (χ, η † ). Thus, the Yukawa interaction becomes After the dark Higgs gets a VEV, the Yukawa interaction generates a fermion mass splitting. The fermion mass matrix is The mass eigenstates are then given by χ 1 = (η + χ)/ √ 2 and χ 2 = i(η − χ)/ √ 2, with masses m 1,2 = m M ± m D . In the mass basis, the dark Yukawa interaction terms can be written as and the interaction with the dark photon is then given by The interaction with the SM is thus off-diagonal, and the direct detection constraint is significantly relaxed because the χ 1 -χ 2 mass splitting means the elastic scattering cross section is suppressed at one-loop (and the one-loop contribution is expected to be small as previously discussed).
In this model, a DM bound state can be produced at the LHC through the process shown in Fig. 1, supported by the exchange of either a dark Higgs or a dark photon. We will focus on the case where the dark Higgs is light and supports the bound state, while the dark photon is heavier and is the principal mediator to the SM, in order to ensure a one-loop suppression in the direct detection cross section while maintaining a large coupling between the quarks and the mediator to the SM and a sizable branching ratio of the bound state to leptons. The dark Higgs is assumed to have some small mixing with the SM Higgs that allows it to decay.
Because of the symmetry breaking pattern, there are only three independent parameters among {m χ , m V , α D , y D }.
The mass hierarchy required above can be achieved by choosing m D m M , so that m D is the small mass splitting, and m 1,2 m M . The spectrum of particles in this model is shown in Fig. 4.
A large value of α D leads to a Landau pole in a broken U(1) theory at a scale above m V [76]. However, since we are mainly interested in the phenomenology of bound states below the scale m V , we assume that a UV completion of the model will avoid the Landau pole. We will later discuss another model with a non-abelian gauge group in the dark sector which will avoid the need for a UV completion.
If the dark bound state, B, is produced from SM initial states, it must be produced from a Z or V exchange. Since the couplings of these gauge bosons to the dark Majorana fermions are off-diagonal, the resulting bound state must be composed of a χ 1 and a χ 2 particle, and for an s-wave state, it must have spin-1. Moreover, since h D only couples χ 1 to χ 1 and χ 2 to χ 2 , decays into final states containing only h D are forbidden, and if m V > m B , the only available decay modes for B are through V back into the SM particles.

D. SU(3)D Triple Higgs Model
We now consider a dark sector model based on a completely broken SU(3) D gauge theory, where all of the phenomenologically desirable properties of the dark sector emerge from the breaking pattern of the gauge symmetry. This model has some similarities with the non-abelian DM models of [77], featuring small mass splittings among the components of the DM multiplet that suppress the direct detection cross section. Because the mediator supporting the bound state is a vector in this model as opposed to a scalar in the pseudo-Dirac model above, the indirect detection constraints of the two models turn out to be quite different.
A completely broken SU(3) gauge group was chosen to allow for a sufficiently large gauge coupling, which is favorable for the production of bound states that are supported by gauge bosons. 3 A broken U(1) theory, such as the one found in the pseudo-Dirac model, with a coupling strength α D 0.5 at momentum scales above the gauge boson mass quickly runs into a Landau pole. Thus a broken U(1) theory with a large coupling constant is likely to have emerged from a larger, nonabelian gauge group in the first place [76]. We choose an SU(3) gauge group rather than SU(2), because for a fermion in the fundamental of a completely broken SU(2) theory with an off-diagonal coupling to the SM, the gauge boson corresponding to the diagonal generator produces a repulsive potential between the two components of the fermion, making it difficult for a phenomenologically viable bound state to exist without introducing additional light mediators.
As in the previous model, the coupling between the dark sector and the SM is mediated by the mixing of the 3 DM models with an unbroken gauge group are constrained by the fact that dark matter is effectively collisionless in galactic dynamics [78].
dark and SM gauge bosons; in this non-abelian case, the mixing operator is non-renormalizable. Bound states in this model are supported by the exchange of one of the SU(3) D gluons, which acquires a relatively small mass during the symmetry breaking. The dark sector contains a triplet of Dirac fermions χ = (χ 1 , χ 2 , χ 3 ) charged under SU(3) D , with a Dirac mass, m χ . After symmetry breaking, the components acquire a small mass splitting, so that m χ1 < m χ2 = m χ3 , with χ 1 and χ 2 ultimately forming an s-wave, spin-1 bound state, B, which can be produced at colliders. χ 1 , being the lightest fermion in this theory, serves as our DM candidate.
The SU(3) D breaking occurs via three Higgs-like fields: two scalars in the adjoint representation of SU(3) D , H 1 and H 2 , and another scalar in the fundamental, H 8 . The dark sector Lagrangian is given by where V µν is the SU(3) D field strength of the dark gluons, a = 1, · · · , 8 is an SU(3) D index and τ a ≡ λ a /2 with λ a being the Gell-Mann matrices.
µ H c i for the two adjoint Higgs fields. The structure of the Gell-Mann matrices is such that V 1 , V 2 couple χ 1 to χ 2 , V 4 , V 5 couple χ 1 to χ 3 , and V 6 , V 7 couple χ 2 to χ 3 ; V 3 couples diagonally to χ 1 and χ 2 , while V 8 couples diagonally to all three fermions; the interaction vertices are shown in Appendix C. The scalar potential V (H 1 , H 2 , H 8 ) can be chosen to satisfy the symmetry breaking pattern that we will describe below.
We impose a Z 2 symmetry at the renormalizable level under which H a 1,2 → −H a 1,2 . This forbids any marginal interaction terms between the Higgs sector and the fermion sector, including a Yukawa interaction term. Therefore, we can treat both sectors as decoupled to first order. However, the following dimension-5 operator is allowed: so that after H 8 acquires a suitable VEV, a mass splitting occurs among the components of χ. Finally, we introduce the following operators that encapsulate the mixing of the dark sector with the SM: (25) Notice that the first term introduces a small breaking of the Z 2 symmetry. This term can originate from a dimension-6 operator that respects this discrete symmetry, such as φH a 1 V a µν B µν , with φ being a scalar field that is odd under Z 2 , which acquires a VEV as well. The details of the origin of this operator are unimportant, as we will focus instead on the phenomenology resulting from the kinetic mixing. 4 At the point of symmetry breaking, H 1 and H 2 acquire a VEV v 1 and v 2 in the 1-and 2-component respectively, and H 8 acquires a VEV given by H 8 = v 8 (cos θ, 0, sin θ), with v 8 m χ v 1 , v 2 and some arbitrary angle θ. This symmetry breaking pattern can be accomplished by choosing an appropriate Higgs potential, which we discuss in detail in Appendix C. The VEV of H 1 in the first term of L mix leads to the conventional kinetic mixing term discussed above, with ≡ 2v 1 /Λ 1 , and V 1 as the mediator to the SM. The second term in L mix guarantees the prompt decay of the other dark gluons through small mixings into the SM: details are discussed further in Appendix C. The choice of H 8 gives a small mass splitting to the Dirac fields in χ, leading to the following fermion masses: We will always neglect the mass splitting when not considering its role in suppressing the direct detection of DM, so that m χ1 m χ2 = m χ3 m χ . The lightest fermion χ 1 is the DM candidate and it is stable; the other particles in the theory decay promptly. More details are provided in Appendix C.
Finally, the dark gluons remain approximately diagonal after the symmetry breaking, with squared masses (up to order g 2 D v 2 8 g 2 D v 2 1,2 ) given by: m 1 also receives O( 2 ) corrections from the kinetic mixing with Z, which we will neglect as was explained above. Thus, the dark gluon masses satisfy the hierarchy and V 8 serves as a good candidate for a bound state mediator. Fig. 5 illustrates the spectrum of particles in this model. 4 One can in principle include the interaction term H a 2 V a µν B µν , but this term does not affect the main features of this model. With the symmetry breaking pattern discussed later, the gauge bosons V 1 and V 2 couple to the same dark fermions, χ 1 and χ 2 . We will leave this term out from the Lagrangian for simplicity. As in the Majorana case, if the dark bound state B arises from SM processes, then it must be produced from the mediator V 1 ; the resulting bound state must be χ 1 χ 2 or its antiparticle equivalent B. Again, since the mediator is spin-1, s-wave bound states must be in the spin-triplet configuration.
In the mass basis, the interaction term responsible for the production is (all fields now denote their mass eigenstate) with r ≡ m Z /m 1 . With m 8 < m χ and the other gluons being significantly more massive than m χ , B's are mediated by V 8 through the interaction terms which leads to an attractive potential between the constituents of B. The coupling between V 8 and the fermions in B is therefore α B = α D /12. The mass hierarchy of this model forbids decays into any of the dark gluons V a for a = 1, · · · , 7. Furthermore, the decay of B into any number of V 8 is forbidden by the conservation of the SU(3) D color charge in the unbroken SU(3): V 8 only couples χ 1 to χ 1 and likewise χ 2 to χ 2 , and cannot carry away the net color charge of B.

IV. EXPERIMENTAL CONSTRAINTS ON DARK SECTOR MODELS
In this section, we will first discuss in Sec. IV A the range of viable model parameters in each of the dark sector models detailed above. We will then study the phenomenology of each of these models at the LHC in Sec. IV B, and their cosmology and indirect detection signatures in Sec. IV C.

A. Viable Model Parameters
In both models, the bound state χ 1 χ 2 (or its antiparticle equivalent, if applicable) is formed from a stable dark matter candidate χ 1 and an unstable fermion χ 2 . In order for the decay of χ 2 to not dilute the production of the bound state, we must ensure that the decay width of χ 2 is much smaller than the decay width of B. In both models, χ 2 decays through an off-shell SM mediator to χ 1 and two SM particles. In the pseudo-Dirac model, this three-body decay width is parametrically Γ χ2 ∼ 2 g 2 D g 2 SM (∆m) 5 /m 4 V , where g SM is a coupling constant to the SM which depends on the actual SM particle considered, and ∆m = m χ2 − m χ1 , which we always take to be small. On the other hand, the bound state decay width is Γ B ∼ 2 g 2 D g 2 SM m 2 χ |ψ(0)| 2 /m 4 V . The relative ratio of these widths is therefore where α B = y 2 D /4π ∼ O(0.1 − 1) for situations where LHC production of bound states is important. An identical relationship holds for the triple Higgs model, with α B = α D /12.
In both of the models we have presented, the interaction between the dark sector and the SM is controlled by a single vector boson: V in the pseudo-Dirac model of Sec. III C and V 1 in the triple Higgs model of Sec. III D. The mixing of the SM and dark sectors shifts the Z mass, and is thus constrained by EW precision tests (EWPT). In particular, the ρ parameter is shifted by an amount [71] where m V is the mass of the SM mediator in either model, and t W is the tangent of the weak mixing angle. The global fit for the central value of ρ to EWPT data is ρ 0 = 1.00037 ± 0.00023 [79]. Constraints are set by requiring that any choice of leads to a minimum value of m V such that ∆ρ is consistent with the 2σ limit for the value of ρ 0 . Next, in order for a bound state to be possible, the constraint given in Eq. (1) must be satisfied. This condition can be satisfied by ensuring that the mass of the particle supporting the bound state is sufficiently small. For the pseudo-Dirac model, this means choosing a sufficiently small dark Higgs mass such that y 2 D m χ > 21.1m h D , and for the triple Higgs model, ensuring that α D m χ > 20.16m 8 .
Finally, to avoid direct detection constraints, the mass splitting must exceed the typical kinetic energy of DM in the solar circle. Taking the velocity dispersion of DM to be v ∼ 10 −3 , this means that the mass splitting has to exceed approximately 10 −6 m χ . A small mass splitting, albeit large enough to be consistent with this lower bound, can be achieved by picking suitable values for the Dirac bare mass m D in the pseudo-Dirac model and Λ m in the triple Higgs model.
In both theories, there are two parameters (m D and m h D for the pseudo-Dirac model, m 8 and Λ m for the triple Higgs model) that can be set to naturally satisfy both the criterion for bound states and avoid direct detection constraints, while having little impact on the LHC phenomenology. However, these parameters can have some influence on the relic abundance of DM in these theories, as well as on indirect detection bounds. This will be discussed after the next section.

B. LHC Phenomenology
We now turn our attention to the production and detection of bound states at the LHC for both theories. In the perturbative picture, bound states B are produced by quark anti-quark parton interactions through an schannel V and Z (mass-eigenstate) boson, with the only available decay mode of B being an off-shell V or Z back into SM particles, leading to resonance signatures. The more accurate procedure of taking into account the mixing of V and B yields a qualitatively similar result; we use the full mixing calculation in all of the plots shown, but focus our qualitative discussion primarily on the perturbative picture. 5 The mono-X + MET search can be effective in setting constraints on these dark models, particularly in the range of parameter space where 2m χ < m V , the region of interest for both dark sector models. To study the constraints that mono-jet + MET searches can place on our models, we use FeynRules [80] and MadGraph [81] to obtain the MET distribution for a wide range of m χ and m V . The distribution is then compared to the observed 95 % confidence upper limit on the number of mono-jet + MET events in 10 inclusive MET bins obtained by AT-LAS with 36.1 fb −1 of data [82]. Any value of m χ and m V with a MET distribution that has more events in any inclusive bin than the 95% upper limit is deemed to be ruled out by the experiment.
Next, we recast bounds from a search for resonance in dilepton events in 36.1 fb −1 of 13 TeV ATLAS data [83] to set constraints on the production of B. In the models considered here, B decays entirely into SM particles with a significant branching ratio to pairs of leptons, making the dilepton resonance search a particularly powerful probe. This search constrains the production cross section times branching ratio of a Z boson assuming some minimal vector couplings to the SM fermions, which allows us to directly interpret these constraints as a limit on the production of cross section times branching ratio of the bound state B.
These searches are also sensitive to the resonant production of the vector mediator V itself, which tends to be significantly more constraining than mono-jet + MET searches when the coupling of the mediator to SM quarks are comparable to the coupling to DM. However, in portal models like the ones we are considering, the mixing into the SM is small while the coupling to DM α D can be large. In the range of parameter space where the mediator mass m V 2m χ , V overwhelmingly decays into χ 1 χ 2 or χ 2 χ 1 , which correspond to final states with MET and are vetoed in dilepton resonance searches to suppress W and Z backgrounds [83]. The search for B, however, faces no such limitation in this region of parameter space.
The production cross section of B (and equivalently of V ) can be computed from Eq. (5), assuming the narrow width approximation. In the perturbative picture, B decays through an -suppressed coupling to the Z, or through V , which has an -suppressed coupling to both J µ EM and J µ Z . The resulting expression for the bound state width to quarks is where α is the EM fine structure constant, Q is the electric charge of the quark, g V and g A are the vector and axial couplings of q to the Z respectively, given by g Note that we have assumed throughout that the bound state is well-approximated by non-relativistic quantum mechanical results, which is a valid assumption so long as the binding energy of B is far less than m χ . For this bound state, we thus require where α B = y 2 D /4π for the pseudo-Dirac case, and α B = α D /12 for the triple Higgs model.
As we argued earlier, the production cross section of B crucially depends on the total width of V ; this means that the total width of V should be included in the computation of the width shown in Eq. (33). Importantly, the width of V should be evaluated at s = m 2 B , since B lies below the χχ open production threshold [79]. The perturbative partial widths of B as well as V into all possible SM final states are shown in Appendix D.
In the mixing picture, the partial widths of V calculated here correspond to Γ V,0 . We take Γ B,0 = 0, since Γ B,0 = Γ χ2 Γ B , as shown in Eq. (31). In the pseudo-Dirac model, there is only one bound state, and the mixing calculation proceeds in the same fashion as described in Sec. II B. The sum of the perturbative partial widths of B, calculated in Appendix D, is numerically a good approximation to the width after mixing, Γ B . For the triple Higgs model, there are two bound states, B and B, and so all three states need to be simultaneously diagonalized. However, B and B maximally mix to form two CP eigenstates, Since V 1 is a CP-even state, it does not mix with the CP-odd combination B − , and the diagonalization is performed over V 1 and the CP-even B + ; the CP-odd state B − does not interact with the SM. In both models, the unmixed mass matrix given in Eq. (3), with the mixing parameter f given by [22,23,26] where N f = 1 for the pseudo-Dirac model, and N f = 1/ √ 2 for the triple Higgs model, which accounts for the differences in coupling and bound-state mixing.
In both models, B cannot decay into final states that only contain the mediator which supports the bound state: this is because both the dark Higgs in the pseudo-Dirac model and V 8 in the triple Higgs model have couplings with the DM fermion number that conserves the number of each of χ 1 and χ 2 .
Decays of B into dark sector final states become possible once m B m V in the pseudo-Dirac model, or m B m 1 in the triple Higgs model: the final states are V h D and V 1 V 8 V 8 respectively. Because of the large coupling between the DM fermions and the mediators, these dark sector decays are the main decay modes of B, rendering the dilepton resonance search for B ineffective. These dark sector final states all mix with the SM, and can in principle lead to multilepton signatures at the LHC. Earlier studies have exploited this signature to look for bound states [11,12], but we do not explore this possibility here for two reasons. First, once m B becomes significantly greater than the SM mediator mass, the resonant enhancement derived in Eq. (8) becomes ineffective, and the cross section for producing B drops quickly away from m B ∼ m V or m 1 . Furthermore, the branching ratio of these mediators to leptons is small, since they kinetically mix through the U(1) Y and decay predominantly into quarks. Second, a direct search for the V or V 1 resonance is significantly more constraining, since the mediator is lighter than the bound state, and there is one fewer factor of the branching ratio to leptons to contend with. In both models, the mediator dilepton resonance search rules out all of the parameter space for m B > m 1 or m V once the coupling to the SM is sufficiently large.
At tree level, we are therefore only interested in the decay modes of B and V into the SM: both particles can decay into a pair of SM fermions, as well as W + W − and Zh where h is the SM Higgs, through the mixing of V with Z/γ. Neither particle can decay into ZZ or γγ final states, since these processes are forbidden by charge conjugation symmetry.
The sensitivity of the dilepton resonance search depends strongly on the width of the resonance, and the 13 TeV ATLAS limits with 36.1 fb −1 of data as a function of the ratio of the width of the resonance to its mass Γ/m are presented in [83]. The total widths of both states are fully taken into account when computing the limits of the search, and the search is assumed to be completely ineffective once Γ/m > 0.32.
The resulting 95% confidence limits from mono-jet + MET, dilepton resonance and EWPT are shown in Fig. 6 in the m χ − m V plane for the pseudo-Dirac model and in Fig. 7 in the m χ − α D plane for the triple Higgs model.
The dilepton resonance search results presented in both figures are searches for the lighter resonance state in the mixing picture; the search switches from V to B along the line m V,0 = m B,0 , where the lighter resonance changes rapidly from being mostly V 0 to mostly B 0 as one moves from below to above this line. 6 As we argued earlier, since the mass eigenstates are always wellseparated and the lighter resonance is always narrow, we can simply assume that the total cross section is given by a Breit-Wigner profile with a width given by either the V (m V,0 < m B,0 ) or the B (m B,0 > m V,0 ) and neglect interference effects. In both cases, the search for the B resonance when m V,0 > m B,0 extends the reach of experimental constraints significantly into this region of parameter space, as compared to what we might expect from just the vector resonance search and the mono-jet + MET search combined.

C. Freezeout and Indirect Detection
We now turn our attention to the freezeout process for the DM in each model, as well as constraints derived from indirect detection experiments. Let us focus on the annihilation channels that do not suffer a suppression by , in order to be as model-independent as possible.
In the pseudo-Dirac model, the potential kinematically available final states (at late times) are h D h D and B h D , with the latter channel corresponding to radiative formation of a bound state, B (which may be spin-1 or spin-0). The V h D final state is forbidden, since V couples χ 1 to χ 2 , and h D couples χ 1 to χ 1 . In the triple Higgs model, if all the gauge bosons and Higgses except V 8 are heavier than the DM, the only open final states are V 8 V 8 and the radiative bound state formation. Note that in the limit where the DM is slow-moving, radiative bound state formation requires not merely that the mediator be light compared to α B m χ , as required for a bound state, but that it satisfy the stronger condition that the mediator mass is smaller than the binding energy, m Y α 2 B m χ /4. Thus, this process can be forbidden by increasing the mediator mass, and indeed we will see that indirect detection limits are much easier to satisfy in regions of parameter space where α B m χ m Y α 2 B m χ /4. In this regime, the DM annihilation products will thus be determined by the decays of the bound state mediator.
During freezeout, the partner particles χ 2 (in the pseudo-Dirac model) and χ 2 , χ 3 (in the triple Higgs model) are also present, and their annihilation and coannihilation channels may also relevant.

Pseudo-Dirac Model
If m V > 2m χ and the bound-state mediator is too heavy for radiative bound state formation, then the only annihilation channel not suppressed by or kinematically forbidden is annihilation to h D h D . Both χ 1 χ 1 and χ 2 χ 2 pairs can annihilate in this fashion, but there is no tree-level coannihilation; χ 1 χ 2 → h D h D does not occur for the same reason that the χ 1 χ 2 bound state does not decay into the dark sector. The cross section for χ i χ i → h D h D in the limit of low DM velocity, before accounting for the Sommerfeld enhancement, is given by: where x h ≡ m h D /m χ . We will assume that during freezeout the mass splitting between χ 2 and χ 1 , set by m D , is small compared to the freezeout temperature; for O(TeV) DM this corresponds to requiring a mass splitting at the GeV scale or below, which is not in tension with the requirement that the mass splitting be large enough to prevent elastic scattering in the present-day halo (where typical kinetic energies for a TeV DM particle are of order 1 MeV or less). In this case, the abundances of χ 1 and χ 2 remain equal during freezeout, as their equilibrium abundances are equal and their annihilation channels are identical. Consequently, each of χ 1 and χ 2 must constitute half the DM abundance, with   the χ 2 subsequently decaying to χ 1 (this occurs through emission of an off-shell V ).

������-����� ��% ���������� ������
Since p-wave processes can dominate during freezeout, to compute the full rate we will need the Sommerfeld enhancement factor for higher-l processes. The Sommerfeld enhancement for multipole l due to a Yukawa potential can be numerically approximated by [44]: We determine the relic density by numerically solving the Boltzmann equation (following the method of [85]) for the χ 1 state and then doubling the result to account for the contribution from χ 2 . We integrate the Sommerfeld-enhanced velocity-dependent cross section over the thermal velocity distribution (assuming a Maxwell-Boltzmann distribution) for the DM at each timestep. As discussed earlier, we neglect radiative bound state formation during freezeout. We define overclosure to occur when Ω χ h 2 > 0.1228, corresponding to the 2σ upper limit (0.1198 + 2 × 0.0015) from Ref. [84].
For m h D smaller than the binding energy, we also account for radiative formation of χ 1 χ 1 bound states (followed by decay into SM particles). To estimate the bound state formation rate via light scalar emission at low velocities, we add to this rate the analytic low-velocity estimate of [51] for the cross section for capture into the ground state (which dominates the overall capture rate), Note that Ref. [51] derives this expression from the Hulthén potential, so in Eq. (40) we have replaced m φ /(α D m D ) in their result with the parameter * φ ; the Hulthén potential with this rescaled mass parameter gives a better approximation to the Yukawa potential [44]. Furthermore, we have included an extra factor of 1/2 to account for the fact that our annihilating particles are identical fermions, and thus only spin-singlet configurations contribute to this s-wave process (yielding a factor of 1/4), but the overall cross section is increased by a factor of 2, as discussed in Ref. [19].

������-����� ��% ���������� ������
The experimental sensitivity to this cross section will depend on the final state to which the h D particles eventually decay, which in turn depends on m h D and whether h D mixes with the SM-Higgs. However, in general hadronic decays will dominate the signal (due to the larger number of hadronic degrees of freedom), and the photon spectra from decays to different quark species are rather similar, as they arise largely from the decays of neutral pions produced in hadronic showers [90]. Thus, we can estimate the sensitivity of indirect detection by examining the constraints set by assuming a bb final state.
In the left panel of Fig. 8 we show limits on the annihilation cross section to bb for Majorana DM from the Fermi [46] and H.E.S.S. [47] gamma-ray telescopes, and sample results for the predicted annihilation rate from our two models. The H.E.S.S. limit, which dominates for DM masses above 1 TeV, is based on a study of the region within 300 pc of the Galactic Center, and assumes an Einasto density profile for the dark matter; if the Milky Way possesses a large core, these limits might be substantially weakened. The Fermi limits are based on a study of Milky Way dwarf spheroidal galaxies. The intermediate step of light mediator production will further broaden the photon spectrum, but Ref. [91] demonstrated that the effect on the constraints is modest for hadronic final states where the spectrum is already quite broad. Thus to obtain an estimate of the constraints, we simply adopt the cross section limits for annihilation to bb. We compare the maximum allowed cross section σv rel max to the predicted cross section scaled by the fraction of DM in the χ 1 state, σv rel (Ω χ1 h 2 /0.1198) 2 ; examples for the pseudo-Dirac model with y D = 2.5 and the triple-Higgs model with α D = 3.0 are shown in Fig. 8, both for m V (m 8 ) = 50 GeV.
In the left panel of Fig. 9, we plot the regions in m χ − m h D plane where bound states exist, the universe is not overclosed, and indirect limits are not violated. We see that there are almost no indirect constraints for DM masses below a few TeV and m h D larger than the binding energy (when m h D is below the binding energy, there remain allowed regions, but they must be chosen to avoid resonant Sommerfeld enhancement). We also plot the regions allowed by indirect detection bounds if a non-thermal history is assumed to ensure that χ 1 constitutes 100% of the DM, with Ω χ h 2 = 0.1198. In this case, the indirect constraints are much more stringent, but the bulk of the region where m h D exceeds the binding energy is still unconstrained.

Triple Higgs Cross Section Limits
FIG. 8: Comparison of predicted DM annihilation rates (including Sommerfeld enhancement and radiative bound state formation) to constraints on thebb channel from Fermi observations of dwarf galaxies [46] and H.E.S.S. observations of the Galactic center region [47]. The red solid line indicates the predicted cross section, rescaled by the fraction of DM squared, for thermally produced DM. For the total DM abundance we take Ωχh 2 = 0.1198 [84]. The red dashed line shows the predicted cross section only, corresponding to the assumption that the annihilating species constitutes 100% of the DM. The region to the right of the vertical purple line is ruled out by overproducing the DM abundance. The left panel shows the result for the pseudo-Dirac model with yD = 2.5; the right panel shows the result for the triple-Higgs model with αD = 3.0.

Triple-Higgs Model
If the vector bosons other than V 8 are all at a heavy mass scale, then the dominant DM annihilation process (not involving bound states) both during freezeout and in the present day is tree-level annihilation to two V 8 bosons. This channel is available forχ i χ i , where i = 1, 2, 3; if σ i denotes the cross section forχ i χ i → V 8 V 8 , then we have: This channel furthermore experiences an attractive swave Sommerfeld enhancement, which for purposes of this estimate we approximate using Eq. (12). Potential exchanges of V 8 bosons, which have large rates compared to processes involving the heavier gauge bosons, do not couple theχ i χ i andχ j χ j states for i = j. Likewise, there is no (tree level) coannihilation to the V 8 V 8 final state. Thus, we can treat the χ i species as evolving independently from each other, annihilating only with their own antiparticles, each experiencing its own long-range attractive Yukawa potential due to V 8 exchange. The effective couplings are α D /12 for χ 1 and χ 2 and α D /3 for χ 3 .
However, one important question is whether the different χ i fields truly evolve independently, and in particular, whether decays and scatterings that interconvert between the χ i states are rapid enough to keep the various state populations coupled during freezeout. An example process isχ 1 χ 1 ↔χ 3 χ 3 scattering via t-channel V 4 or V 5 exchange (see Appendix C). As all such processes involve the heavier gauge bosons, they are slow compared to annihilation into a V 8 V 8 final state near the time of freezeout. For the χ 1 χ 1 → χ 3 χ 3 process, the cross section for this scattering process is approx- D /m 2 χ , the rate of processes that scatter one type of fermion into another is suppressed by a factor of ∼ m 4 χ /M 4 . Thus, the process χ 1χ1 → χ 3χ3 freezes out before the χ 1χ1 → V 8 V 8 and is therefore not relevant to determining the relic abundance. This estimate assumes that all of the gluons other than V 8 are more massive than the DM; if this assumption breaks down, the three dark-matter-like populations will no longer evolve independently, and freezeout will be modified.
Under this assumption, we solve separate Boltzmann equations for each of the χ i species (accompanied by their antiparticles), and require that the resulting mass density 2(m χ1 n χ1 + m χ2 n χ2 + m χ3 n χ3 ) matches the cosmological density of dark matter. The masses of the three states are assumed to be equal, with mass splittings small compared to the temperature at freezeout. The greater annihilation rate ofχ 3 χ 3 causes its abundance to be depleted faster thanχ 1 χ 1 andχ 2 χ 2 .

������ ����� �������� ��������� ������
To estimate the late-time indirect detection limits, we proceed as for the pseudo-Dirac case above, and present our results in the right-hand panels of Figs. 8-9. The allowed cross section for DM annihilation is doubled as the DM χ 1 is a Dirac fermion in this case. Since there is an unsuppressed s-wave annihilation channel, there are useful constraints from indirect detection even when radiative bound-state formation is kinematically forbidden. To estimate the contribution from bound-state formation, we numerically calculate the cross section for capture into the ground state by dipole photon emission, and also add the contribution from an s-wave initial state transitioning into the first excited state by emission of a dipole photon. The former process dominates when the mediator mass can be neglected [19], but is suppressed in the very-low-velocity regime as it corresponds to a pwave initial state [92]; thus we add the latter process to properly include the leading contribution at very small velocities. We follow the numerical method described in Ref. [19].
Note that in this case, the scenario where χ 1 constitutes 100% of the DM at late times is essentially completely excluded by indirect detection, for α D = 3.0 and m χ below 10 TeV; such a scenario requires a non-thermal origin for the DM, as the annihilation cross section is well above the thermal value and would deplete the DM density efficiently during freezeout. If non-thermal processes produce more DM at late times, then the large bare annihilation cross section and accompanying Sommerfeld enhancement (and possibly radiative bound-state formation) gives rise to very strong indirect detection signals, as shown in the right panel of Fig. 8.
Both the overdepletion of the DM density and the large direct detection signals may be avoided if the Diracfermion DM possesses some tiny asymmetry, similar to the baryon asymmetry of the SM. The large annihilation cross sections found in these models can readily deplete the DM abundance to the point where the asymmetry sets the residual relic density, and then the indirect-detection signals are suppressed by the absence of the symmetric component. Note that if no such asymmetry is present, indirect detection limits may also pose challenges for sub-GeV DM and mediators as studied by Ref. [11]; thermal relic DM can be quite generically excluded for sub-GeV mediators and sub-TeV DM [93].
This behavior does not occur for the pseudo-Dirac model because the principal annihilation channel is pwave suppressed; this both makes it possible for TeVscale χ 1 particles to constitute 100% of the DM with a thermal history, and ensures that large regions of parameter space remain that are not excluded by current indirect detection bounds (although bound state formation can provide indirect detection signals, as in Ref. [51]).

V. DISCUSSION
The resonant production of dark sector bound states at the LHC can be an important complementary search channel to the missing energy and mediator resonance searches. Unlike a mediator resonance search, a bound state resonance search directly probes the properties of the DM, and can be more effective when the mediator decays primarily to invisible DM particles. In addition, a bound state resonance search can be more sensitive than a missing energy search strategy at high DM masses.
We have studied the general features of models that can be probed by bound state resonance searches at the LHC while remaining consistent with other powerful experimental constraints. These models generally require a sufficiently strong coupling to a light mediator that can support the bound state, and a heavy mediator that couples the dark sector to the SM. We also carefully take into account the mixing between the heavy mediator and the dark matter bound state. Bound state decays into the light mediator should also be suppressed to allow for a significant partial width into SM particles.
These requirements must be reconciled with constraints from both direct and indirect detection experiments. Spin-independent direct detection cross sections can be suppressed by having only loop-level interactions between the dark sector and nucleons, which can result from an off-diagonal coupling between the SM and two DM states with a mass splitting between them. Constraints from gamma-ray experiments and overclosure must be carefully considered, taking into account Sommerfeld enhancement due to the presence of a light mediator and radiative bound state formation both during freeze-out and in the present day.
The SU(2) L minimal DM models possess many of the properties that we have discussed above, but pure wino or higgsino DM chargino bound states have a production cross section that lies well below the sensitivity of dilepton resonance searches, although DM particles in larger representations of SU(2) L with a large electric charge forming a deeply-bound electromagnetic bound state can potentially be discovered.
We propose two dark sector models with kinetic mixing into the SM that contain bound states that can be probed effectively through bound state resonance searches at the LHC, while remaining consistent with direct and indirect detection constraints. The pseudo-Dirac model contains two Weyl fermions with a small mass splitting between them, capable of forming bound states through a light Higgs mediator, while the triple Higgs model is an SU(3) gauge theory with a single Dirac fermion in the fundamental representation, with all of the properties required for a viable model being generated by symmetry breaking of the gauge group. We study the LHC phenomenology of these models and find that bound states searches are complementary to both missing energy and vector mediator resonance searches, and are particularly powerful at high DM masses. A simple rescaling of our constraints indicates that future 27 TeV or 100 TeV pp colliders could potentially probe DM bound states with masses in the O(10) TeV range.
We find that these models naturally avoid overclosure of the universe, and broad swaths of parameter space exist where they also evade limits from indirect detection searches under the assumption of a thermal history, despite the presence of the bound state implying modelindependent large enhancements to the low-velocity annihilation rate. The indirect limits are most easily satisfied when radiative capture to the bound state in the local DM halo is kinematically forbidden, because the mass of the mediator supporting the bound state exceeds the binding energy.
If the bound-state-forming species is required to constitute 100% of the DM, through a non-thermal history, but is symmetric in the present day, then the indirect searches are sufficient to rule out almost all of the parameter space of interest in the LHC bound-state resonance search for the triple Higgs model; the pseudo-Dirac model evades this fate through a late-time velocity suppression of its annihilation rate. Where a DM species has a greaterthan-thermal annihilation cross section but still constitutes 100% of the observed DM density, a viable model that can be first detected by resonance searches at the LHC should possess some suppression to the annihilation cross section at late times, due e.g. to a dominant p-wave annihilation channel or a small primordial asymmetry.
To summarize, dark sectors with bound states can be probed at the LHC through resonance decays to SM particles. Models with multiple force carriers and DM-like states, where the DM scatters inelastically off SM quarks at tree-level, can naturally give rise to a sufficiently large production cross section while evading direct detection constraints. The presence of a light mediator, needed to support the bound state, modifies freezeout and leads to stringent indirect detection limits; however, these constraints leave a wide region of parameter space open, while suggesting a preferred mass spectrum where the mediator mass exceeds the binding energy. DM models with bound states possess a rich phenomenology, allowing complementary probes from many different experimental directions. the draft. We particularly thank Yevgeny Kats for pointing out the importance of mixing between bound states and mediators. Feynman diagrams in this paper were generated with TikZ-Feynman [94], and Feynman amplitude calculations performed with FeynCalc [95,96] where V µν is the dark sector gauge field strength. To diagonalize this, we define the new field B µ = B µ + V µ , and thus B µν = B µν + V ab. µν , where ab. indicates the abelian part of the field strength tensor. All terms of O( 2 ) are neglected here. Then The last term is an additional interaction that is unimportant in the diagonalization. In terms of the new field, we have the SM fields where we have defined Z µ = c W W 3 µ − s W B µ and A µ = s W W 3 µ + c W B µ . In terms of these new fields, the mass terms for Z and V are: Diagonalizing this and defining r ≡ m Z /m V (we neglect the O( 2 ) shift in the masses), we get the following mass eigenstates to first order in (marked by tildes): with A being the massless mode. With this, we see that the original SM fields become and This is the result shown in Eq. (17).  In the triple Higgs model of Sec. III D, H 1 and H 2 acquire a VEV v 1 and v 2 in the 1-and 2-directions respectively, at a scale above the bare Dirac mass m χ . This breaks the gauge symmetry down to a residual U(1). Subsequently, this remaining symmetry group is broken at a lower scale than m χ by H 8 , which acquires a VEV given by with v 8 m χ . This symmetry breaking pattern can be achieved with the following Higgs potential, which obeys the Z 2 symmetry mentioned above: The last two terms forbid a VEV in the second component of H 8 , which would make H † 8 τ a H 8 non-zero for a = 1. This symmetry breaking pattern produces the mass hierarchy for the dark gluons and the mass splitting between the Dirac fields in χ, which were discussed in the main text.
After symmetry breaking, the kinetic mixing terms with the SM become where we have defined 8 /2 ≡ v 2 8 /Λ 2 8 , taking 8 1. The first term represents the kinetic mixing between dark sector and the SM discussed in the main text, while the remaining mixing terms are highly suppressed but nonzero for generic values of θ: their existence guarantees that the gluons V 3 , V 4 and V 8 decay to SM particles over cosmological timescales.
With the m 8 < 2m χ < m 1,··· ,7 mass hierarchy, the only stable dark sector particle is χ 1 , since we can assign a conserved dark baryon number to the fermions and χ 1 is the lightest dark fermion. The heavy gluons V i for i = 1, · · · , 7 can decay into a pair of dark fermions since their masses exceed 2m χ . Decays into a pair of dark fermions for V i with i = 2, 5, 6 and 7 occur promptly for an O(1) coupling g D . For V 8 , which mixes directly into the SM and can decay into a pair of SM fermions, its lifetime is approximately where ∆m χ ≡ v 2 8 /2Λ m is the mass splitting between χ 2 , χ 3 and χ 1 . For any reasonable choice of parameters considered in this paper, it is clear that V 8 , χ 2 and χ 3 are all unstable on cosmological timescales, and that the only dark matter component in the present day is χ 1 . Table I shows the perturbative partial decay widths of the dark sector bound state B in both models into SM final states, through mixing with the SM mediator V . The bound state wavefunction is given in Eq. (34). In addition, Table II shows the perturbative partial decay widths of the SM mediator V into all possible final states.

Appendix E: Bound State Formation via Initial/Final State Radiation
In addition to the resonant formation process that has been our main focus in the body of this work, bound states can also form in conjunction with radiation of other particles in the initial or final state. This process is very important in the context of electron accelerators where the center-of-mass energy of the colliding particles is fixed and does not overlap the bound state resonance (as discussed e.g. in Refs. [11,97]), and so the resonant signal is absent.
This process could also be critical if the decays of spin-0 bound states were much more observable than those of spin-1 bound states, and the mediator with the SM were a vector (or if the spin-1 states were more observable and the mediator were a scalar); emission of additional particles would then allow the production of the rarer but more observable bound state. However, in the examples we have studied, the latter situation does not hold; indeed, the spin-0 s-wave bound states can generally decay into light mediators and are thus likely to be more difficult to detect than their spin-1 counterparts.
Since initial and final state radiation inevitably involves extra powers of the coupling relative to the resonant case, we expect processes of this type to be suppressed relative to the resonant production. However, one might wonder whether threshold enhancements to the production and interaction cross section for unbound but slow-moving DM particles, in the presence of a light mediator, could modify this conclusion and lead to a large contribution from the threshold region.
Note that this is a very different parameter regime to that considered for muonium production in Ref. [97] and for light darkonium production in Ref. [11], where the beam energy is presumed to be large relative to the resonance energy, and the extra particle(s) emitted as initial/final state radiation carry away much of the beam energy; it is more similar to the situation in indirect detection, where slow-moving DM particles may emit a light particle and radiatively capture into a bound state (see e.g. [41,50]). The rate for such radiative bound-state formation scales as 1/v close to threshold, for a massless mediator. However, we will show that in the case where the particles are produced near threshold and then form a bound state, this 1/v scaling is canceled out by the small phase space for the particle production near threshold.
Similar contributions to bound-state formation from soft gluon emission have been studied in the context of quarkonium formation using non-relativistic effective field theory techniques [98,99]. In that case, p-wave color-singlet quarkonia can be formed either directly or through an intermediate s-wave color-octet pair of heavy quarks; relatedly, the s-wave quarkonium state |QQ can be described as having a small O(v 2 ) admixture of a Fock state |QQg containing an additional soft gluon. This approach suggests that the admixture term can be neglected to leading order when dealing with s-wave bound states, and should not experience large enhancements near threshold.
To see explicitly how this works in our case, note that we can write the matrix element for production of the bound state (plus a light mediator with momentum l), via an intermediate state of two near-threshold (i.e. highly non-relativistic) but unbound DM particles, as: HereM(i → a 1 , a 2 ) is the hard matrix element describing production of two free DM particles with momenta a 1 , a 2 from the initial state i, and likewiseM( b 1 , b 2 → q 1 q 2 l) is the hard matrix element describing the radiation of a light mediator with momentum l from the DM-DM state with particle momenta b 1 , b 2 , to produce final-state DM particles with momenta q 1 , q 2 . The wavefunctions convert the plane-wave states to the full intermediate and final states accounting for potential effects. p 1 and p 2 act as labels on the intermediate state with momentumspace wavefunctionψ p1,p2 , describing the momenta of the constituent particles at large separation.ψ B denotes the momentum-space wavefunction of the bound state (which in principle is labeled by the quantum numbers n, l, m; we suppress these indices). m χ is the DM mass and M ≈ 2m χ is the bound-state mass. In the non-relativistic limit where the potential is neglected, the leading-order matrix element for light vector boson radiation from one of a pair of heavy fermions (with gauge coupling g B and fermion masses m 1 , m 2 ) is given by: Inserting this expression into Eq. (E1), setting the masses of the two heavy fermions equal, m 1 = m 2 = m χ , working in relative momentum coordinates, and choosing the center-of-mass frame, we obtain: × d 3 q (2π) 3ψ * B ( q) q ψ p1,p2 ( q + l/2) +ψ p1,p2 ( q − l/2) .
The integral over d 3 q on the last line also appears in the matrix element for radiative bound state formation, and has been previously computed in the non-relativistic limit for massless vector mediators [19,100]. In the nearthreshold regime, l α 2 m χ (as the binding energy must provide the necessary energy to radiate the mediator), and the l-dependence of the integral can be neglected; in this case, the integral simply scales as 1/ √ p, where p = ( p 1 − p 2 )/2. (This factor, when squared, is responsible for the 1/v scaling of the radiative bound state formation cross section.) If we further suppose that the hard matrix element for production of the intermediate state from the initial state is independent of the final-state relative momentum a, i.e. we can write iM(i → a 1 , a 2 ) = iM(i → DM, DM) then the integral over d 3 a simplifies to give iM(i → DM, DM)ψ * p1,p2 (0), where ψ denotes the position-space wavefunction. The wavefunction at the origin in a Coulomb-like potential scales as α B m χ /p (e.g. [100]), which yields the usual Sommerfeld enhancement when squared.
Putting these pieces together and performing the phase-space integral over d 3 p 1 , writing E p1 = E p2 = m 2 χ + | p| 2 since we are working in the COM frame, we find that (keeping only scaling relationships, dropping order-1 factors): iM(i → f ) ∼ g B √ m χ iM(i → DM, DM) Note that as mentioned previously, the phase-space integral over the intermediate-state momentum d 3 p has canceled out the 1/p scaling from the wavefunctions.
Thus we see that the bound-state production cross section through this channel should scale as |M(i → DM, DM)| 2 α 2 B , multiplied by a 2-body phase space factor. Since the momenta in the final state are small, of order l ∼ α 2 B m χ , the overall scaling of the cross section with the couplings is α 4 B × |M(i → DM, DM)| 2 . By comparison, the resonant production cross section scales as |M(i → DM, DM)| 2 α 3 B , where the α B dependence arises from the B wavefunction. Thus the rate to produce an extra light mediator by emission from a near-threshold intermediate state, in conjunction with the bound state formation, is suppressed by one power of α B overall. This is the same suppression one would naively expect for emission of a hard photon from the initial or final state, with no small phase-space factors or threshold enhancements. We self-consistently neglect all such diagrams in the body of this work.
Here we have neglected the mediator mass m Y in estimating the scalings; in particular, the intermediate-state position-space wavefunction may be steeply peaked near the origin for special values of m Y , corresponding to the presence of near-zero-energy bound states (e.g. [40]). However, it seems likely that any apparent enhancement from this behavior can be reinterpreted as resonant capture into a near-zero-energy bound state, which is already accounted for in our formalism. We leave a more detailed study of the resonant regime to future work.