Broadening Dark Matter Searches at the LHC: Mono-X versus Darkonium Channels

Current searches for dark matter at the LHC focus on mono-X signatures: the production of dark matter in association with a Standard Model (SM) particle. The simplest benchmark introduces a massive spin-1 mediator, the $Z^\prime$ boson, between the dark matter $\chi$ and the SM. Limits derived from mono-X channels are most effective when the mediator can decay into two on-shell dark matter particles: $M_{Z'}\gtrsim 2M_\chi$. We broaden the experimental reach into the complementary region, where the $Z^\prime$ mediator is much lighter than the dark matter. In this scenario the $Z^\prime$ mediates an effective long-range force between the dark matter, thereby facilitating the formation of darkonium bound states, as is common in many dark sector models. The darkonium becomes active when $M_{\chi}>M_{Z'}/\alpha_{\rm eff}$, where $\alpha_{\rm eff}$ is the effective fine-structure constant in the dark sector. Moreover, the darkonium could decay back into SM quarks, without producing missing transverse momentum in the detector. Considering multijet final states, we reinterpret existing searches to constrain the simple $Z^\prime$ benchmark beyond the region probed by mono-X searches. Assuming a baryonic $Z^\prime$ mediator and a Dirac dark matter, direct detection bounds can be loosened by giving a small Majorana mass to the dark matter. We also consider the interplay between mono-X and darkonium channels at future high energy colliders, which is at the frontier of probing the model parameter space.


Introduction
One important mission of the Large Hadron Collider (LHC) and future high energy colliders is to probe the nature of dark matter. If the dark matter particle has a coupling to the standard model sector, it could be produced at the LHC, usually in pairs if the dark matter is stabilized by a (possibly new) symmetry. The dark matter particles are expected to escape the detector like neutrinos. They can lead to events with large missing transverse momenta, if another visible object (e.g., an energetic jet) is produced at the same time. The monojet process has been widely studied at the Tevatron, LHC and future colliders [1][2][3]. The same idea has been extend to other standard model particles being produced in together with dark matter leading to the so-called mono-X searches [4].
Going beyond mono-X, another important aspect of dark matter at colliders is the production of dark bound states. Bound states made of dark matter and its anti-particle (darkonium) exist generically in dark sector models with a dark force carrier whose coupling to the dark matter is strong enough. They are the analog of the positronium or heavy quarkonium states in the real world, which have played an instrumental role in our understanding of the SM. It is conceivable that similar phenomena would occur in a dark sector containing the dark matter [5][6][7][8][9][10][11][12][13]. The signatures of darkonium have been studied at both lepton and hadron colliders in several models [14][15][16][17]. The bound state formation channel is also an ideal place for probing the self-interactions of dark matter in the laboratories [15].
In this work, we investigate the complementarity between the mono-X and and darkonium channels in the LHC search for dark matter. Our study is based on a simple renormalizable model where a Dirac fermionic dark matter χ is charged under the gauged baryon number symmetry. The new Z boson mediates the interaction between dark matter and quarks. This simple model is widely adopted as the benchmark for LHC monojet analyses [18,19]. So far, experimental limits have been derived in the region of parameter space with M Z 2M χ , where the Z can decay into two on-shell dark matter particles. Outside of this region the production rate of a pair of dark matter particles through the off-shell Z is too small and the mono-X searches become ineffective. It is possible to directly search for the production of Z which subsequently decays back into the SM quarks, resulting in multijet final states [20][21][22][23][24]. However, the resonance search in the multijet final states quickly loses its constraining power for a Z at or below the weak scale, due to the overwhelming QCD background. Therefore, there is presently no experimental search that is sensitive to the Z benchmark when the Z is light.
In this work we would like to point out that, in the commonly adopted benchmark for mono-X searches, the Z boson could mediate a long-range dark force between dark matter particles, when its mass is light and coupling to dark matter strong. Then the χχ darkonium bound states could exist in nature and be produced at a high energy collider. Once produced, the χχ inside the darkonium will eventually find each other and annihilate, causing the latter to be unstable and decay back to SM quarks. The novelty here is that, although the dark matter particle is produced at the collider, there is are missing energy/momentum in the final state! In this case, the darkonium would appear as a resonance in multijet final states and its production can be constrained in these searches. In turn, an experimental limit on the production rate of darkonium can be translated into constraints on the mass and couplings of the dark force carrier: the Z . In the end we find the darkonium signals are most active when the Z -quark coupling is weak and the Z -dark-matter coupling is strong. When darkonium exists, it offers a new handle to explore the nature of dark matter at colliders, and can be highly complementary to the mono-X channel as well as the direct searches for the Z boson.
This paper is organized as the following. In section 2, we describe the simple benchmark model and discuss the necessary condition for the darkonium bound states to exist, which includes requiring the Z to be lighter than the dark matter, precisely the region where the mono-X search in ineffective. We give a brief summary on the existing searches for a light baryonic Z boson. In section 3, we calculate the darkonium production cross section and the possible decay channels. We explore the feasibility of using the di-jet channel to search for the darkonium states appearing as new resonances. We derive the existing LHC limit as well as the projections at the future high-energy high-luminosity LHC, as well as a possible 100 TeV pp colliders. These results are compared with the reach of the monojet channel. We highlight the complementarity of mono-X versus darkonium searches, both of which are needed to effectively cover each other's blind spot. In section 4, we discuss the implications from other areas of dark matter searches, including direct and indirect detections, as well as the its production mechanism in the early universe. We identify the parameter space where high-energy colliders are at the frontier of searching for dark matter in this model. Then we conclude in section 5.
While this paper was being prepared, a related work [25] appeared which explored dark matter bound state signals in several non-minimal dark sectors with quite sizable dark couplings. However, the simple benchmark model discussed in this work was not covered.

The Benchmark Model
In mono-X searches the commonly adopted simplified model includes a massive spin-1 boson, the Z , mediating the production of the dark matter particle, which is assumed to be a vector-like pair of fermions (χ,χ). The leading low-energy effective Lagrangian takes the form The Z is assumed to have a universal coupling g q to SM quarks and, to be general, we allow for both the vector and axial-vector current couplings with the dark matter. The axial coupling then implies the Z current is anomalous, which can be remedied by postulating spectator fermions to restore the gauge invariance associated with the Z . We further assume these spectator fermions to be much heavier than the weak scale. Since we will not consider loop-induced processes involving the Z boson in this work, the anomalous Z current (or equivalently the anomaly-cancelling spectator fermions) plays no role in our study [26].
A concrete example of a Z boson is to gauge the baryon number symmetry U (1) B in the SM, which is anomalous with respect to the electroweak gauge groups. The ultraviolet complete models of gauged baryon number have been discussed in [27][28][29]. In Eq. (2.1) we have also extended the minimal gauged baryon number model by introducing an additional dark matter field χ that is charged under the U (1) B . The presence of the axial coupling g χ implies that χ L and χ R must carry different charges under the U (1) B . As a result the dark matter mass M χ is not U (1) B invariant and must be generated via the Yukawa coupling of χ to the vacuum expectation value (vev) of a U (1) B Higgs field. The same vev also contributes to the mass of the Z . In the appendix A, we present a simple model for this. Requiring the Yukawa coupling to satisfy the unitarity constraint results in an upper bound on g χ (see Eq. (A.6) and also [30]), Therefore, if M Z M χ , which is the region of interest in this work, g χ need to be small. For a dark matter interacting with the quark through a Z mediator, constraints from dark matter direct detection are quite stringent for a light Z [31][32][33][34]. Such constraints could be relaxed by introducing, in addition to the Dirac mass term for χ, a Majorana mass in Eq. (2.1) [35][36][37][38], We assume δ is small enough compared to the Dirac mass M χ so that our discussions on collider phenomenology in section 3 remain valid at the zeroth order in the small δ expansion, which allows us to treat χ as a Dirac fermion in collider studies. On the other hand, δ must be large enough to evade the direct detection constraints. A quantitative estimate of δ satisfying both considerations will be presented in section 4. A non-vanishing δ will have implications in cosmology and indirect detection of dark matter, which will also be explored in section 4.

The Formation of Darkonium
One important aspect of dark matter we want to explore is bound state physics. The Z exchange yields a Yukawa potential between χ andχ. With a light enough Z and large enough couplings g χ , g χ , bound states made of χ andχ could form. Because of the fermionic nature of χ, there are two darkonium ground states, one with total spin S = 0 and the other with S = 1, which we denote by η D and Υ D , respectively. We will focus on ground states in this work. The vector coupling of the Z with the dark matter yields an attractive Yukawa potential while the sign of the potential from the axial coupling depends on the total spin [39]. We can define the effective fine-structure constant of the Z -mediated long range interaction between the dark matter particles as The potential is attractive for S = 1 and repulsive for S = 0. The Z boson plays two roles in this model. It is not only the mediator between the dark matter and the SM, but also the dark force carrier responsible for self-interactions of the dark matter.
At the LHC, the darkonium can be created via an off-shell Z boson, much like the production of J/Ψ particle through an off-shell photon in QCD. Therefore, the spin-1 darkonium Υ D can be singly produced on resonance, while the spin-0 darkonium η D has to be produced in association with another Z . * In what follows we will focus on the spin-1 darkonium Υ D , in which case Then the condition for the ground state to exist is [40], The mass of Υ D is given by 2M χ minus the ground state binding energy, BE. In the Coulomb limit (M Z → 0), BE = α 2 eff µ/2, where µ is the reduced mass of the system For general nonzero M Z , the binding energy can be solved numerically [40]. A useful analytic approximation can be obtained using the Hulthén potential to mimic the Yukawa potential [41]. In this case, where a 0 ≡ 1/(α eff µ). One could also derive the bound state wavefuction at the origin, which is [41], The single production of Υ D at the LHC could be described, effectively, by a kinetic mixing with the Z boson, which takes the form [15] Through this kinetic mixing Υ D could couple to SM quarks. Any non-zero axial current coupling will introduce further kinematic mixings of both the Z and Υ D (known as the 1 −− ground state) with the 1 ++ state, an excited bound state made of χχ. However, these mixings are suppressed compared to Eq. (2.10) by additional powers of α eff . We will therefore truncate the spectrum and only consider the ground state for the rest of this paper.

Current Constraints on the Z Mediator
Experimentally, a vector boson Z that couples to SM quarks could be produced at hadron colliders such as the Tevatron and LHC. If the Z is lighter than twice of the dark matter mass, it can only decay back to a SM quark and antiquark. Existing dijet resonance searches cover the Z mass window from 50 GeV to multiple TeV scales. Below we list several limits from the recent analysis on the Z -quark-qntiquark coupling from dijet searches (see also [24] These limits directly apply to the Z in our model when it predominantly decays into qq † . The future running of LHC and the dijet searches could further improve the coverage of Z mass from 50 GeV up to a few TeV, leaving the region of light Z below 50 GeV as a blind splot. Through the quark loops, the Z boson mixes with the SM Z-boson. As a result, when the masses of the two are close enough, there are useful limits from the hadronic Z-boson width measurement at the LEP [42,43]. For even lighter Z , below a few GeV, there are also constraints on its mixing with the heavy quarkonium states like the Υ and J/Ψ, as well as the rare decay of meson states into the Z . For a recent study, see [44]. It is also worth noting that for very light bayonic Z , the heavy anomalon fields can have strong non-decoupling effects on flavor-changing neutral current processes [45,46]. However, our study here will mainly focus on the region M Z > 10 GeV, thus these non-decoupling effects can be evaded.
In the next section, we will show that the same search results could be reinterpreted as constraints on the production of the darkonium Υ D , leading to new limits on the dark matter simplified model that are complementary to the mono-X searches.

Darkonium Versus Mono-X
In this section, we will explore the interplay between darkonium and mono-X channels in searches for dark matter whose interaction is mediated by a U (1) B baryonic vector boson. They turn out to be highly complementary to each other in probing the model parameter space. Moreover, the dijet Z search at LHC seems to have a blind spot for light Z below 50 GeV. As explained above, a Z is light enough could facilitate the existence of darkonium bound states. Search for the formation of such new states at the LHC could in turn constrain the light Z as a dark force and help covering the above blind spot. These important features are summarized in Fig. 1. We will go through the details of this plot for the rest of this section. Generically, the mono-X searches are most sensitive to the region M Z 2M χ , while the darkonium searches mainly probe the region M Z α eff M χ .
Hereafter, we will choose the follow benchmark values for the model parameters, Before moving on, we comment on how our results change when the benchmark values of α χ and g q vary from the choice in Eq. (3.1). First, from Eq. (3.2) (see below), the dark matter bound state production cross section is proportional to α 4 χ . The monojet cross section (dominate by ISR jet radiation) is proportional to α χ . Their limits will get substantially weaker for smaller values of α χ , especially for the bound state channel. We also restrict ourselves to the region with α χ < 1 so that the χ particles in the bound state are still non-relativistic and we could reliably do perturbative calculations in the small α χ expansion. Second, we choose a relatively smaller value of g q than commonly used in the previous monojet analysis (where g q 0.2 is used). This is mainly driven by the increasingly stronger bound from the Z search in the dijet channel. With a coupling g q 0.2, most the region in Fig. 1 with M Z > 50 GeV is already excluded by the current LHC data.
We will scan the rest of parameter space and present our results in the M Z versus M χ plane. We find that, with the current LHC data (13 TeV, ∼ 36 fb −1 ), we are not yet able to derive a competitive limit in the parameter space of interest. However, future experiments such as the upcoming high luminosity running of LHC at 14 TeV (expected luminosity up to ∼ 3 ab −1 ), the high-energy highluminosity LHC running at 27 TeV (expected luminosity up to ∼ 15 ab −1 [47]), and a possible 100 TeV hadron collider [48], will enable us to derive very useful limit in the parameter space where the dark matter bound states could be produced. In Fig. 1, we show the region of parameter space that could be probed by these future experimental programs.

Mono-X Searches
The simple model in Eq. (2.1) has served as the benchmark model for many mono-X searches for dark matter at the LHC, where the dark matter particles χ andχ are produced in together with a SM particles. In particular, the "monojet" final states are characterized with very large transverse missing energy (MET) plus one or more jets. A representative Feynman diagram is shown in Fig. 2 (left). The qg initiated process is dominant because of the large gluon parton distribution function (PDF) at small x. In the parameter space where M Z > 2M χ , the Z boson could be produced on-shell in association with one or more jets, qg → q + Z , followed by the decay Z → χχ resulting in MET.
Recent monojet analyses of the model by ATLAS and CMS collaborations can be found in Refs. [18,19]. With g q = 0.1 and α D = 0.5, constraints derived from rescaling the results in Refs. [18,19] do not yet place useful limits on the parameter space shown in Fig. 1. However, according to our estimates, this will change with the future running of the LHC at higher luminosities (300 fb −1 and 3 ab −1 ). The reaches are shown by the cyan dot-dashed and dotted curves in Fig. 1. The region to the left of these curves could be covered. With a large enough luminosity, the monojet constraint extends slightly into the region where M Z < 2M χ , where the monojet production is a 2 → 3 process, qg → qχχ, with the Z being off-shell. Nevertheless, the cross section decreases rapidly with increasing dark matter mass, because the radiated jet needs to have a large transverse momentum, of order a few hundred GeV, to satisfy the experimental trigger. This feature limits the ability of using monojet channel to probe the parameter space deep in the M Z < 2M χ region.
Instead of initial state jet radiation, one may also consider final state radiation of the Z boson, qq → χχZ . In the M Z < 2M χ region, the Z can only decay back to qq, which appear as two jets. For a sufficiently light and boosted Z , the two jets will be collimated with each other and may appear as a single jet in the detector. In this case one could apply the monojet analysis to this channel.  However, the final state radiation process must be initiated by qq initial states, see Fig. 2 (right), and the cross section is suppressed by the anti-quark PDF over the gluon PDF compared to initial state radiation case. We include this channel in our analyses and find the modification to the total monojet cross section to be small (less than 10%). It is possible to study this channel further by exploring the ⌥ D < l a t e x i t s h a 1 _ b a s e 6 4 = " 3 y a Z 7 4 q 8 o z 3 6 L 1 4 r 9 P W n D e b 2 Y V v 8 N 4 + A e g n k J s = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " 3 y a Z 7 4 q 8 o z 3 6 L 1 4 r 9 P W n D e b 2 Y V v 8 N 4 + A e g n k J s = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " 3 y a Z 7 4 q 8 o z 3 6 L 1 4 r 9 P W n D e b 2 Y V v 8 N 4 + A e g n k J s = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " 3 y a Z 7 4 q 8 o z 3 6 L 1 4 r 9 P W n D e b 2 Y V v 8 N 4 + A e g n k J s = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " 3 y a Z 7 4 q 8 possible jet substructure [49], as well as displaced vertex [50] signatures.

Darkonium Searches
The limitation of mono-X searches outside the M Z < 2M χ region strongly motivates us to consider additional possible dark matter production channels at the LHC, in particular, bound states of χ,χ. These states are unstable and will decay promptly back (the decay rates are given by Eq. (3.3)) to SM quarks, appearing as dijet (or multi-jet) resonances, which lead to very different collider signatures from monojet. We want to emphasize again that the darkonium search here is different from direct searches of Z as dijet resonances discussed in Section 2.2. Here the Z plays the role of a dark force for the darkonium to exist. The darkonium production cross section is proportional to its wavefunction at the origin thus depends on the Z mass and couplings. Constraining the formation of such darkonium state allows us to indirectly constrain the dark force. We discuss these in detail in this subsection.

Υ D production
As discussed in section 2.1, we will focus on the spin-1 darkonium state Υ D at the LHC. It is mainly produced via qq fusion and the Feynman diagram is shown in Fig. 3 (left). The production cross section at a proton-proton collider takes the form Here we calculated the cross section using the NNPFD [51] with the PDF set NNPDF30 lo as 0118 nf 6. After the production, Υ D will decay into two (or more) jets as will be discussed in the section 3.2.2. We will use the dijet resonance search data to set limits and estimate future reach at the LHC and higher energy colliders.

Υ D decay
After production, there are three ways for the darkonium Υ D to decay: 1) to qq via an off-shell Z ; 2) to two Z bosons; 3) to three Z bosons. The partial decay rates are where N f is the number quark flavors that Υ D can decay into, and Ψ(0) is given in Eq. (2.9). The calculation of non-relativistic bound state decay is reviewed in [52]. The first decay channel is simply the inverse process of the Feynman diagram in Fig. 3 (left). The second decay channel Υ D → Z Z is possible only in the presence of nonzero g χ coupling, which violates the charge-conjugation (C) parity. The Feynman diagram for this process is shown in Fig. 3 (right). The effective operator responsible for this decay channel is [53]Ô For the Υ D → 3Z decay rate, we work in the limit that g χ g χ and M Z M χ . This allows us to derive an analytic expression for the decay rate, in analogy to that of Υ → 3γ decay in the SM [54]. Using the value of g χ from Eq. (3.1), we find that Γ Υ D →2Z Γ Υ D →3Z , i.e., the three-Z decay is always subdominant.

Dijet resonance search for Υ D
With the above production and decay channels, we are now ready to quantify the experimental constraints for Υ D by recasting dijet resonance searches. To date, the ATLAS and CMS collaborations have published several results on the dijet resonance search [20][21][22][23], which covers the resonance mass from ∼ 50 GeV to multiple TeV scales. These searches assume the heavy resonance to have 100% decay branching ratio into qq.
However, in our model, Υ D , as the heavy resonance, has more than one decay channels. In order to properly interpret the LHC limits from dijet resonance searches, we need to simulate the selection efficiency of each possible decay channel of Υ D in Eq. (3.3). To this end we first create a FeynRules [55] model containing both the Z boson and the spin-1 darkonium Υ D . In the model file we include the kinetic mixing in Eq. (2.10) responsible for the production of Υ D , as well as the effective coupling in Eq. (3.4) that mediates the Υ D decay. Then we use MadGraph 5 [56] to generate the Υ D production and decay to jets at pp colliders, and run PYTHIA 8 [57] and DELPHES 3 [58] for hadronization and detector simulations. We follow the dijet event selection cuts described in [20][21][22][23] to derive the efficiency factor, Eff i , for each Υ D decay channel. In Fig. 5, we show the efficiency factors for the Υ D → qq and Υ D → 2Z → 2(qq) channels to pass the event selection cuts in each mass window, which are called Eff 1 and Eff 2 , respectively. We simulate the production of Υ D at the LHC and take into account of its boost on event-by-event basis. A lighter Υ D is typically born with a higher boost, thus when it decays the opening angle of final states tends to be smaller, leading to a lower efficiency factor. Such an effect is shown in the upper left panel of Fig. 5.
It is worthwhile remarking on the dijet efficiency factor for Υ D → 2Z → 2(qq) decay, which is the following. Kinematically, when M Υ D M Z , the Z bosons from the decay of Υ D are boosted.
For an Υ D produced at rest, the two jets from each Z have a maximal opening angle The formation of darkonium requires M Z ≤ 6α eff M χ /π 2 ≈ 0.6 α eff M χ , which leads to (∆R jj ) max 0.6 with the benchmark parameters. If the two jets are within the cone size of θ jj < 0.4, they will be reconstructed as a single jet typically. There is an order 1 chance for this to occur. This estimate is confirmed by the plots in the second row of Fig. 5. We find it convenient to define the effective coupling between Υ D and SM quarks where i goes through all the possible Υ D decay channels labelled in Eq. (3.3). The sub-label "1" stands for the Υ D → qq decay channel. The upper limit on the effective coupling g Υ can be directly read from the existing LHC limits on elementary Z -quark-antiquark coupling obtained in [20][21][22][23], for four mass windows (which is called g q there). Because g Υ D is a function of the all model parameters in Eq. (3.1), an upper limit on g Υ D will translate into a contour in the parameter space in Fig. 1.
We find that the current LHC data are not yet able to provide a competitive constraint in the plot. However, the further running of high energy high luminosity LHC (at 27 TeV), as well as the possible 100 TeV collider will do. To estimate the future reaches, we first scale the number of events with the increasing integrated luminosities, by a factor where L now are given in [20][21][22][23]. We then calculate the enhancement factors in the production cross sections for both the signal, and the background, and consider √ s CM = 14, 27, 100 TeV as the future collider energies. The Υ D production cross section is given by Eq. (3.2). The QCD background cross section for dijet production at parton level goes as, ∼ŝ −1 . Note that the dijet search is a bump hunt. In practice, we focus on a narrow dijet invariant mass windowŝ ∼ M 2 Υ D . As a result, the proton-proton level cross sections are proportional to the following quantities (the parton luminosity defined in [59]), respectively where τ = M 2 Υ D /s CM . We calculate the rescaling factors R sig √ sCM and R bkg √ sCM using the NNPDF.
Therefore, the future upper bound on g Υ D is expected to get stronger by a factor of The future collider reaches are shown in Fig. 1 for three of the mass windows (blue, red, green curves, with texts next to them denoting the corresponding future collider energy and luminosity). The regions below these curves could potentially be covered.

Impact of the Majorana mass term on collider phenomenology
So far, our discussions of collider phenomenology are based on the effective Lagrangian Eq. (2.1) but with the Majorana mass term for χ defined in Eq. (2.3) set to zero. Here we clarify the impact of a nonzero δ on the dark matter spectrum and the bound state physics LHC. In the presence of both M χ and δ, the mass terms for χ can be written as Here we assume δ is a real parameter. The two Majorana fermion mass eigenstates and the corresponding eigenvalues are In terms of χ 1,2 fields, their interaction terms involving the Z boson now take the form (3.14) In the parameter space of interest to bound state physics, M Z M χ , the constraint Eq. (2.2) indicates that the diagonal axial-current interactions are suppressed, g χ g χ . The vector-current interactions are dominant and they must be off-diagonal with respect to χ 1,2 . In this case, bound states made of a χ 1 and a χ 2 particle can still form [38]. The long range force due to Z exchange alternates the two states along each fermion line (see Fig. 6).  One can still write down a Schrödinger equation describing such a bound state, with the reduced mass now defined as Using this reduced mass instead of that in Eq. (2.7), one can repeat the discussions in section 2.1 to find the spectrum and wavefunctions. In the small δ/BE expansion, where the binding energy BE is defined in Eq. (2.8), the two results must agree at the leading order. In Fig. 7 we show in green the region of parameter space with δ > BE. Outside of the green we could have δ BE ∼ α 2 χ M χ . For example, we show in green dashed line in Fig. 7 where δ/BE = 1/10. We expect the main results on bound state collider phenomenology, derived based on pure-Dirac fermion assumption in the previous subsections, remain unaltered.
This said, in the presence of nonzero δ, the χ 2 particle becomes unstable. For δ > M Z , the following decay could occur, χ 2 → χ 1 Z , whose decay rate is For Λ QCD δ < M Z , the decay of χ 2 has to occur through off-shell Z , χ 2 → χ 1 qq. In the case δ M χ , M Z , the decay rate takes the approximate form For δ Λ QCD , the final state qq will turn into meson states. Isospin singlet vector mesons can directly mix with the baryonic Z boson. The decay rate for χ 2 → χ 1 ω is where f ω 70 MeV is the decay constant of the ω meson, ω|ūγ µ u +dγ µ d|0 = √ 2f ω m ω ε µ ω . For δ < m ω , χ 2 could decay into χ 1 plus pions via off-shell ω; and for δ < 2m π , χ 2 has to decay into χ 1 plus e + e − (or µ + µ − ) through the (loop generated) kinetic mixing between Z and the photon. In practice, we require χ 2 must not decay within the time scale of the bound state formation, which is equivalent to requiring Γ χ2 to be smaller than the bound state binding energy. For this reason, in Fig. 7, we also shade out the region with Γ 2 > BE in blue color.

Direct, Indirect Detections and Early Universe
In this section, we discuss the implication of dark matter direct and indirect detection constraints on the model parameter space which was explored in the previous section, using the same set of benchmark parameters given in Eq. (3.1). We also address the possible (thermal) origin of our dark matter relic abundance from the early universe.

Direct detection
We first consider dark matter direct detection, in the presence of a nonzero δ parameter. In this case, the dark mater splits into two Majorana mass eigenstates, χ 1 and χ 2 . Without loss of generality, we assume χ 1 is the lighter one and exist in nature as the dark matter. χ 2 is an unstable partner state. With the Z -quark-antiquark coupling in Eq. (2.1) and the Z -χ 1 -χ 1,2 couplings in Eq (3.14), there are two types of χ 1 -nucleus scattering processes. One is spin-independent and inelastic, χ 1 + N → χ 2 + N , whose cross section is proportional to the product of couplings, g 2 χ g 2 q . The other process is elastic, χ 1 + N → χ 1 + N , whose cross section is proportional to g 2 χ g 2 q , and depends on the spin of χ 1 . Because the SM quarks still couple to Z coherently via their number density, the spin vector of χ 1 has to be contracted with either its velocity v, or the three momentum transfer q. As a result, the cross section is also velocity dependent and suppressed by the halo velocity squared (v halo ∼ 10 −3 c). This suppression makes the latter cross section safely small in view of the current direct detection limits.
Next, we will examine the inelastic scattering more carefully. The nucleus-level scattering cross section in the small δ limit is where T is the target nucleus, and µ 1T = M χ1 M T /(M χ1 + M T ). We also assume M Z is much larger than the momentum transfer of the scattering. The state-of-art dark matter direct detection limits are obtained by the PandaX-II [60], LUX [61] and XENON1T [62] collaborations, where for dark matter mass of a few hundred GeV, the upper limit on the nucleon-level cross section is σ SI 10 −45 cm 2 . The nucleon level scattering cross section can be calculated as which is much larger than the current upper bounds. ‡ If this was the case, most of the parameter space shown in Fig. 1 would have been ruled out, where we explored the LHC searches for dark matter.
The only way to suppress this cross section and get around the constraint is to turn on δ. The phase space factor in Eq. (4.1) implies a minimal χ 1 velocity for the scattering to occur, v ≥ v min = 2δ/µ 1T [35,36]. The usual assumption is that the halo dark matter velocities satisfy the Maxwell-Boltzmann distribution which is peaked around v peak 270 km/s and has a cutoff at the escape velocity v esc 544 km/s [63]. Therefore, if the minimal velocity v min v peak , the population of χ 1 that could trigger the scattering process is exponentially suppressed, and if v min > v esc the process will be turned off completely.
We have calculated the lower bound on δ numerically so that the direct detection limits are satisfied, which is shown by the orange region in Fig. 7. Approximately, this bound coincides with the kinematic limit, The main message from Fig. 7 is that there exist a large window of δ (the white region) where our collider discussions remain valid and the direct detection constraints are evaded.

Thermal relic abundance
Next, we discuss the dark matter relic abundance in this model. We will make the most modest assumption that the dark matter χ 1 and the SM particles were in thermal equilibrium with each other in the early universe. Its relic abundance is obtained thermally via the freeze out mechanism. There are several ways for χ 1 to annihilate in the early universe. When M χ1 > M Z , two χ 1 particles can annihilate into two Z bosons via a t-(or u-) channel χ 2 exchange. The annihilation cross section is given by When M χ1 < M Z , the above annihilation channel is forbidden, unless one or both of the Z bosons goes off-shell. We take into account another channel where two χ 1 particles annihilate into qq via an s-channel off-shell Z . This is only possible via the diagonal Z coupling in Eq. (3.14) which is an axial-current interaction involving two χ 1 particles. Its cross section is given by, Here, the annihilation cross section is P -wave suppressed. A simple way of understanding the Pwave nature is from parity. The total parity of a fermion-anti-fermion system (applies to two χ 1 particles) is (−1) +1 , where is the orbital angular momentum between the two particles. The axial current (spatial part) has even parity. Therefore we must need = odd for the annihilation amplitude to be non-vanishing. It is worth noting that during the thermal freeze out v rel 6T f /M χ1 and When the temperature is high enough, χ 2 particles also exist in the universe. As a result, Eq. (3.14) permits another annihilation channel, the χ 1 and χ 2 coannihilation. The cross section is , (4.6) which is an S-wave annihilation. For the values of δ (which controls the χ 1 -χ 2 mass splitting) allowed in Fig. 7, we find δ T f . Thus, the relative Boltzmann suppression between χ 1 and χ 2 populations is not significant. In this case, we have (σv) χ1χ2→qq (σv) χ1χ1→qq because we have chosen g χ g χ . There is also another annihilation channel χ 1 χ 2 → Z Z , which involves one diagonal and one off-diagonal coupling from Eq. (3.14), and is proportional to g 2 χ g 2 χ . However, we find this cross section is subdominant to χ 1 χ 1 → Z Z which is an S-wave annihilation contains a g 4 χ term (Eq. (4.4)), again because g χ g χ . In Fig. 8, the blue solid contour shows where χ 1 could obtain the observed dark matter relic abundance [64], by requiring the total annihilation cross section for χ 1 to be equal to σv th 3 × 10 −26 cm 3 /sec. We also draw two contours of constant relic density χ 1 in unit of the observed dark matter relic density (labelled by f 1 = 0.02 and 10 −3 ). We neglect the non-perturbative Sommerfeld corrections to the cross sections which is usually an order one effect for thermal freeze out. The shape the blue contours are similar to those found in [65], although in our model we have also kept g χ non-vanishing thus more annihilation channels have been included. In the light blue shaded regions, the χ 1 annihilation cross section is smaller than σ th , thus the dark matter would be overproduced in a thermal history. Outside the blue shaded regions in Fig. 8, the relic abundance of χ 1 is underproduced. In this case, its relic density in unit of the observed dark matter relic density is § where (σv χ1 ) tot is the sum of the cross sections in Eqs. (4.4)-(4.6). This means that χ 1 can only comprise a fraction of the total dark matter in the universe.

Indirect detection
Next, we examine the dark matter indirect detection constraints, assuming χ 1 comprises (a fraction of) the dark matter candidate. We assume χ 2 do not exist in the universe today. The indirect signals could arise from χ 1 χ 1 annihilation in the galaxy or the early universe. We will take into account of the lower limit of δ derived from direct detection constraints in Eq. (4.3). With a nonzero Majorana mass, one cannot make the assumption that the dark matter relic abundance in the universe is asymmetric and argue away the indirect detection constraints [66][67][68]. The Born level annihilation cross sections included in this calculation are Eqs. (4.4) and (4.5). For the P -wave annihilation, Eq. (4.5), it is worth noting that the dark matter halo velocity is a small number, v rel ∼ 10 −3 , thus this cross section is highly suppressed.
On top of the Born-level cross sections, we also take into account of the possible Sommerfeld effect in the total annihilation rate. This is especially important when the mass of Z is smaller compared the de Broglie wavelength of dark matter. We calculate this non-perturbative factor by numerically solving the Schrdonger equation, following the pioneering works [31,33,[69][70][71][72][73].
One clarification is necessary with a non-zero Majorana mass δ, where the usual Sommerfeld effect derived for pure Dirac fermion case needs to be modified. The key picture is that a long-range Z  Fig. 1 with the same choice of parameters. Assuming a thermal history of the dark matter χ 1 , it could obtain the correct thermal relic density along the blue solid curve. The blue shaded region is because of the overproduction of χ 1 's relic density. Outside the blue region, χ 1 is underproduced and could only account for a fraction of the total dark matter (shown by blue dot-dashed contours). The red spiky regions are ruled out by indirect detection experiments due to the Sommerfeld enhancement, in spite of the small relic density. All the white regions in this plot are still alive. Here, f 1 is the fraction of observed dark matter relic density that is comprised of our dark matter candidate χ 1 , defined in Eq. (4.7). The "monojet" and "darkonium" territories denote the regions of parameter space where the monojet and darkonium resonance search channels at LHC are most powerful, as discussed in section 3 and shown in Fig. 1.
exchange converts the χ 1 χ 1 initial state into χ 2 χ 2 intermediate state. Because the typical potential energy is of order ∼ α 2 χ M χ , the usual Sommerfeld effect only applies for δ α 2 χ M χ . We will assume that this is the case for simplicity. If δ is too large compared to the potential energy, one can no longer cut the ladder diagrams which now becomes genuinely loop suppressed. ¶ The interplay between δ and BE in dark matter self interaction was noticed and explored in detail in [38,74].
With the benchmark parameters given in Eq. (3.1), we can evaluate the effective cross section for χ 1 annihilation today, where S (S ) is the Sommerfeld factor for an S-(P -) wave annihilation process, and the factor f 2 1 takes into account that χ 1 may only comprise a fraction of the observed dark matter relic density in our model, which is derived based on Eq. (4.7). ¶ This also corresponds to the green region in Fig. 7.
An analysis on hidden sector dark matter annihilation has been performed in [75] which takes into account that the annihilation into SM particles (quark and antiquarks here) could occur via multiple steps. We adopt the model independent constraints from there. For the dark matter mass range of interest to this work, the Fermi gamma ray observation from dwarf galaxies [76] gives the strongest upper bound on (σv) eff . In Fig. 8, the red regions show the parameter space which is ruled out by this indirect measurement. The spiky feature is mainly due to the Sommerfeld effect. Clearly, the indirect detection constraint can only exclude very limited regions. In the region M Z < 6α eff M χ /π 2 where dark matter bound states exist, we still need the future running of LHC and higher-energy colliders to effectively probe this region (see also discussions in section 3.2).

Conclusion
The nature of dark matter remains mysterious to us after a tremendous amount of effort in searching for them. This strongly suggests going beyond the existing approaches and cast a wide net. One important aspect is to broaden the mission of existing experiments. In this work, we propose reinterpreting the LHC di-jet (multi-jet) resonance search results to look for darkonium bound states which occur in dark sector models with a light dark force carrier and a sizable dark coupling with dark matter. We focus on a simple model where the dark matter interacts with standard model quarks via the exchange of a vector boson Z , which is the same benchmark model widely employed by mono-X searches at the LHC experiments. In the parameter space where the Z is weakly coupled to quarks but strongly coupled to the dark matter, we show that darkonium channel is most useful and highly complementary to mono-X searches. Both ought to be included and contrasted in the analysis of future results from LHC and higher energy colliders.
We have also considered the dark matter production in the early universe as well as direct and indirect detection constraints. We identify the parameter space where these constraints could be weakened, and the reasons behind. The strong direct detection limits can be evaded by turning on a small Majorana mass for dark matter and split the Dirac fermion into two Majorana particles. As a consequence, this excludes the possibility of accommodating the asymmetric dark matter scenario thus indirect detection must be considered. To derive the population of our dark matter today, we resort to a thermal history and assume it acquires its relic density via the freeze out mechanism. Because the dark gauge coupling of interest to this work is order one (for bound states to exist), for dark matter below a few TeV, it is underproduced and could only comprise a fraction of the observed dark matter relic density. This suppresses the indirect detection limits even in the presence of strong Sommerfeld enhancement effects. We find the above effects occur in a large portion of the model parameter space, where the collider searches is the most powerful in probing dark matter in this model.
Although all our findings are based on a very simple model, it is worth emphasizing that the mono-X versus darkonium complementarity as well as some of the dark matter features derived here are generic and applicable to many extended dark sector models. been authored by Fermi Research Alliance, LLC under Contract No. DE-AC02-07CH11359 with the U.S. Department of Energy, Office of Science, Office of High Energy Physics. The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a non-exclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this manuscript, or allow others to do so, for United States Government purposes.
A Origin of a dark matter axial-current coupling to Z In this appendix, we present a simple model that could generate an axial current dark matter coupling to the Z . Under the gauged U (1) B symmetry, the left-and right-handed components of dark matter χ have different charges. The charge assignment is We assume q L = q R . The Lagrangian that respects U (1) B takes the following form Because χ L and χ R have different U (1) charges, we cannot directly write down a mass term, but instead a Yukawa coupling with φ. We assume the scalar potential V (φ) is such that φ get a non-zero vacuum expectation value, φ = w/ √ 2. This vev breaks the U (1) gauge symmetry giving a mass to V and also give mass to the fermion χ via the Yukawa coupling.
The particle mass spectrum after the symmetry breaking is The gauge coupling between χ and Z can be rewritten as Compared to the definition of parameter we have been using, we have If the charges q L and q R are close to each other, the coupling g χ is suppressed by the difference, so is the mass of the vector boson M Z = 2g χ w. In together with the fermion mass, we also find the relation In general the value of the Yukawa coupling y is bounded from above by perturbative unitarity, roughly y √ 4π. Therefore, we find an upper bound similar to the one given in Eq. (2.2),