Collider Signatures of Coannihilating Dark Matter in Light of the B-Physics Anomalies

Motivated by UV explanations of the $B$-physics anomalies, we study a dark sector containing a Majorana dark matter candidate and a coloured coannihilation partner, connected to the Standard Model predominantly via a $U_1$ vector leptoquark. A TeV scale $U_1$ leptoquark, which couples mostly to third generation fermions, is the only successful single-mediator description of the $B$-physics anomalies. After calculating the dark matter relic surface, we focus on the most promising experimental avenue: LHC searches for the coloured coannihilation partner. We find that the coloured partner hadronizes and forms meson-like bound states leading to resonant signatures at colliders reminiscent of the quarkonia decay modes in the Standard Model. By recasting existing dilepton and monojet searches we exclude coannihilation partner masses less than 280 GeV and 400 GeV, respectively. Since other existing collider searches do not significantly probe the parameter space, we propose a new dedicated search strategy for pair production of the coloured partner decaying into $bb\tau\tau$ final states and dark matter particles. This search is expected to probe the model up to dark matter masses around 600 GeV with current luminosity.


Introduction
The Standard Model (SM) of particle physics is the most accurate description of the fundamental particles and their interactions. There are, however, solid experimental and theoretical reasons to postulate the existence of New Physics (NP), e.g., astrophysical observations have established the presence of dark matter in the universe [1]. Since dark matter cannot be accounted for by the particle content of the SM, its nature remains one of the biggest mysteries in modern physics.
At the high-energy frontier, no definite Beyond the SM (BSM) signals have emerged from the full set of run-II LHC data. However, in the last decade a large number of lowenergy flavour experiments performed at B-factories and LHCb have uncovered indirect hints of lepton flavour universality violation in B-mesons decays in b → s + − and b → c ± ν transitions. Interestingly, these deviations from the SM, known as the B-physics anomalies, seem to point towards a new boson at the TeV scale, which couples predominantly to third generation fermions. The only models that can simultaneously accommodate all B-physics anomalies while satisfying the rich low-energy phenomenology and high-p T constraints feature leptoquarks as the mediators of the NP effects [2][3][4][5][6][7][8][9][10][11][12][13][14][15][16].
In particular, the vector leptoquark U 1 with SM quantum numbers (3, 1, 2/3) has emerged as the only single-mediator solution [17][18][19]. Within the scope of ultraviolet (UV) complete frameworks, this vector leptoquark can arise as a gauge boson of a spontaneously broken gauge group containing SU (3) c colour as a subgroup. The minimal phenomenologically viable group is SU (4) × SU (3) × SU (2) L × U (1) T 3 R ; models based on this group are called '4321 models' [20][21][22][23][24][25][26][27]. These particular UV models predict other heavy gauge bosons, e.g., a Z and a colour octet (coloron), that give rise to a rich set of collider signatures that can currently be tested at the LHC [28,29]. In addition, these models introduce a wide range of new fermionic states, which may also lead to new collider signatures.
Inspired by these UV constructions, we explore the possibility that dark matter is contained within a new fermionic multiplet, which would be intimately connected to QCD through a larger gauge symmetry. 1 In particular, we investigate whether dark matter could be a neutral remnant of a fermionic multiplet of the broken gauge symmetry, much in the same way leptons are colourless remnants of a multiplet unifying quarks and leptons (as in the Pati-Salam model [51]). This would suggest that there may be other dark sector fermionic partners similar in mass to the dark matter particle that carry colour and electric charge. The U 1 leptoquark and other new gauge bosons may then mediate interactions between the dark sector and the SM.
Furthermore, a Majorana mass term generated by higher-dimensional operators induces a pseudo-Dirac scenario [52] where the would-be Dirac dark matter fermion splits into two Majorana fermions, the lighter being a dark matter candidate. We show that the quasi-degeneracy of the dark matter candidate and its coloured partner can be naturally preserved within our general UV framework, and can have a significant impact on the thermal freeze-out of dark matter and collider phenomenology. For example, the relic abundance will be determined by coannihilating effects in the early universe. The coloured coannihilating partner also hadronizes and forms QCD bound states, which may lead to unexpected signatures at the LHC.
The best probe of the dark sector is LHC production of two coannihilation partners. These can either form a bound state with each other, which then decays into dileptons via electroweak interactions, or each partner can decay via a leptoquark, leading to a bbτ − τ + and missing transverse energy signature. Since the b and τ particles will be soft, due to the quasi-degeneracy of the dark sector, existing ATLAS [53,54] and CMS [55,56] searches are not sensitive. We propose a dedicated search in this channel which exploits a novel observable: the ratio between the visible and missing energy of the process. By simulating signal and background events, we estimate the reach at the LHC (with run-II luminosity) of both a cut-based and a multi-variate analysis, finding a significant improvement when using a boosted decision tree (BDT) classifier.
The paper is structured as follows. In section 2 we first motivate the simplified model from a UV perspective, and then present the Lagrangian relevant at the TeV scale. The coannihilation effects are described in section 3 and the parameters required to reproduce the observed dark matter relic abundance are calculated. In section 4 we describe the hadronization of the coannihilating partner and the resulting meson spectroscopy. We then examine the main experimental constraints set by existing collider searches and derive bounds on the relevant parameter space. Finally, in section 5, we present our dedicated search strategy for the pair production of the coannihilating partner.

Theoretical Motivation and the Simplified Model
Although we will study a simplified model, we first motivate Majorana dark matter coannihilating with a slightly heavier colour triplet partner from a UV perspective. Under the SM group SU (3) c × SU (2) L × U (1) Y the dark sector states will have the representations χ ∼ (1, 1, 0) , ψ ∼ (3, 1, 2/3) , (2.1) and couple directly to the vector leptoquark U 1 ∼ (3, 1, 2/3). There are two main ways of accounting for a multi-TeV vector leptoquark: gauge models and strongly interacting models [28]. In this work we will assume a gauge model explanation, where the leptoquark is a massive gauge boson associated with the spontaneous breaking of a gauge symmetry, G NP ⊃ G SM .

General Setup
We assume that in addition to the SM matter content, there is also a heavy vector-like fermion transforming under some representation of G NP that decomposes as under the SM gauge group. Although the simplest example is to take X = (ψ 1 ψ 2 ψ 3 χ) T transforming as the fundamental of SU (4), we emphasise that the setup is more general and eq. (2.2) is strictly all that is required. After the breaking of G NP , ψ becomes coloured under SU (3) c and χ remains colourless and electrically neutral, with quantum numbers given by eq. (2.1). We imagine that X is charged under some stabilising symmetry while the SM particles are not, so the components of X constitute a dark sector and the neutral component leads to a dark matter candidate. The breaking of G NP → G SM around the TeV scale typically requires a set of scalar fields uncharged under the stabilising symmetry that we denote collectively as Ω. Since X is a Dirac field by construction, we introduce non-renormalizable interactions so that the neutral component χ gives rise to Majorana dark matter. 2 This then points towards a higher scale (beyond the TeV) associated with fermion number violation. Before the spontaneous symmetry breaking of G NP , the effective Lagrangian for the dark sector up to dimension d = 5 is where Λ is a large cut-off scale associated with fermion number violation and O (5) n are effective operators with Wilson coefficients c (5) n . Since X is a vector-like fermion, the operators O (5) n fall into two possible categories: (i) operators of the form X F X that conserve fermion number, and (ii) Weinberg type operators of the form X c F X that violate fermion number by two units. Here F and F are d = 2 scalar operators that are bilinears in the fields Ω responsible for spontaneous symmetry breaking. Once these scalar fields develop vacuum expectation values (vevs) the O (5) n terms in eq. (2.3) will lead to both Dirac and Majorana masses for the components of X of the form F /Λ and F /Λ, respectively. These masses will be suppressed because the scale Λ associated with fermion number violation is taken to be much larger than the G NP breaking scale. In fact, we will assume that Λ 2 F , F , m 2 X . In section 2.2 we provide a simple 4321 extension that gives rise to these two types of d = 5 effective operators.
Note that the Lagrangian eq. (2.3) is not completely general since there could also be d = 4 Yukawa terms of the form X ΩX , if the scalar sector Ω contains states with the correct quantum numbers (e.g., the adjoint scalar Ω 15 in the 4321 model discussed in appendix A). Once these fields acquire vevs, X ΩX could lead to a mass splitting of order Ω between the ψ and χ components of the X multiplet, while still preserving colour. We here assume that these type of d = 4 terms are not present or are negligible (e.g., due to small Yukawa couplings) and that the only relevant sources of mass-splitting at tree level between the X components arise from d = 5 operators or higher. Higher order effective operators with d ≥ 6 can also be included in eq. (2.3), but these will have little impact on the dark sector phenomenology and for this reason can be dropped.
After the breaking of G NP , the scenario described above will lead to the following mass terms for the dark sector particles:

4)
2 Although the field content in our setup is similar to that considered in ref. [44], the dark matter candidate in ref. [44] is a Dirac fermion. For Dirac dark matter, the dominant interaction that connects the dark and SM sectors is via a Z gauge boson and the strongest experimental limits come from direct detection experiments. Once appropriate non-renormalizable operators are included, the dark matter candidate becomes a Majorana fermion, the dominant interaction is mediated by a U1 vector leptoquark and direct detection constraints are negligible. Instead, collider searches are the best probes of the motivated parameter space.
where m ψ and m χ are two Dirac masses roughly of the same order, satisfying m ψ = m X and |m ψ − m χ | ∼ F /Λ m X , and m L and m R are two small Majorana masses of order ∼ F /Λ, such that m L,R m X . Note that in general we expect m 2 L,R F as the Majorana mass terms in eq. (2.4) violate the G N P symmetry, so they may only break it softly at scales lower than the breaking scale. The resulting mass eigenstates χ 1,2 are two pseudo-Dirac fermions [52,57] with quasi-degenerate masses m χ 1,2 , given by up to leading order in |m L − m R |/m χ and the Majorana condition χ 1,2 χ c 1,2 is satisfied up to the same order. In the mass basis, the Lagrangian reads At tree level, the spectrum of the dark sector satisfies m χ 1 m χ 2 m ψ with mass splittings of order F /Λ ∼ F /Λ. In what follows we will be interested in the compressed spectrum scenario m χ 1 m ψ m χ 2 , so that the lightest stable state χ 1 is a dark matter candidate. Small deviations from eqs. (2.5) and (2.6) can occur at loop level if additional heavy states couple differently to ψ and χ. For example, the heavy gauge bosons in the 4321 models produce order 10% mass splittings between ψ and χ [44].
In summary, we see that relatively general considerations lead to models with a dark sector containing a coloured partner ψ and pseudo-Dirac pairs χ 1 and χ 2 , where χ 1 , χ 2 and ψ have a compressed spectrum.

UV Realisation: Majorana Dark Matter in a 4321 Model
In order to further motivate the general setup presented above, we present a concrete example of a UV model which gives rise to the dark sector effective Lagrangian eq. (2.3). For this, we extend the matter field content of the standard (or flavoured) 4321 model (see appendix A for more details) with a dark sector containing a Z 2 -odd vector-like fermion X ∼ (4, 1, 1, +1/2) with components X = (ψ, χ) T , and a Z 2 -odd right-handed fermion singlet S R ∼ (1, 1, 1, 0). The role of this exact Z 2 symmetry is to stabilise the dark matter. The Z 2 -even SM fields are exactly as in the 4321 models, as shown in table 3. The scalar sector of 4321 contains the Z 2 -even state Ω 1 ∼ (4, 1, 1, −1/2), which gives rise to Yukawa interactions between the dark sector fields. The Lagrangian of the dark sector reads It is convenient to rotate the singlet fields into the Majorana basis and work with the field N defined by 9) such that N = N c and where the θ phase corresponds to the argument of the Majorana mass, M S = |M S | exp(iθ). We also define the couplings λ R ≡ λ exp(−iθ/2) and λ L ≡ λ exp(−iθ/2). We can now rewrite the terms involving S R in eq. (2.8) as The Majorana mass M N ≡ |M S | is assumed to be much larger than the Dirac mass of X , M N m X , since it is generated at a very high scale where fermion number is violated. This allows us to integrate out the N field, giving rise to the effective Lagrangian given in eq. (2.3) with the following d = 5 operators and corresponding Wilson coefficients: and their Hermitian conjugates. The first two operators O LR conserves fermion number, leading to a Dirac mass term. After Ω 1 obtains a vev, Ω 1 = (0, 0, 0, ω 1 / √ 2) T , the 4321 group breaks down to the SM gauge group and the physical masses in eq. (2.7) are given by Notice that for λ L = λ R this UV model reproduces the desired compressed mass spectrum with the correct ordering m χ 1 < m ψ < m χ 2 . Furthermore, the eigenstates χ 1,2 naturally arise as a pair of pseudo-Dirac fermions since the χ 1 -χ 2 mass splitting terms m L,R = λ 2 L,R ω 2 1 /2M N are suppressed by the seesaw mechanism in the limit M N ω 1 .

Field
Type

The Simplified Model
In our analysis, we will study a simplified model motivated by the above scenario. We will parametrise the mass splittings as focusing on the regime where 0 ≤ ∆ ψ 0.3, as explained in the next section, and where ∆ χ 2 2∆ ψ (motivated by eqs. (2.14) to (2.16)).
The covariant derivative terms in L DS , eq. (2.7), contain couplings to gauge bosons which become massive after G NP breaking. These will act as mediators between the dark sector and the SM. As discussed in the introduction, the B-anomalies suggest the existence of a U 1 vector leptoquark. Once this is included, closure of the algebra means that there must also be a heavy Z [28]. Although there could be further gauge bosons, such as a coloron, we assume that these do not significantly impact the dark sector phenomenology. The complete set of new fields we introduce in our simplified model is shown in table 1. We will see that the relic abundance and LHC phenomenology is driven by χ 1 , ψ and U 1 . In the notation of ref. [58] which classifies co-annihilation models and their LHC signatures, this is model ST2 (note that [58] uses a different hyper-charge convention to the one we employ here). The interaction terms of the new gauge bosons are where the numerical prefactors are motivated by the 4321 model [24]. Note that the Majorana fields χ 1 and χ 2 satisfy χ i γ µ χ j = −χ j γ µ χ i for i, j ∈ {1, 2}, so the vector currents χ i γ µ χ i vanish and the above Lagrangian is Hermitian.
The O(1) coupling g U and the masses m U , m χ 1 , m ψ are the parameters which dominate the determination of the dark matter abundance and collider physics signatures we consider.
For these parameters, the key relation that provides a solution to the B-physics anomalies is [29] g For the β and ζ couplings we work in the flavour basis where the SU (2) L SM fermion doublets are aligned with the down-quark sector, where V is the CKM matrix. In this basis, motivated by the B-anomalies and minimality, we assume the following structure for the β couplings Non-zero values for β ij L are needed to address the B-anomalies, and may arise when connecting this simplified model to UV models addressing the structure of the SM Yukawa couplings [25,27]. While these parameters only have a marginal influence on the dark matter abundance and the collider physics signatures we consider in this work, for definiteness we take β sτ L = −β bµ L ≈ 0.11 and β sµ L = −β dτ L ≈ 0.02 [25]. For the Z couplings, we take ζ 33 q = ζ tt u = ζ bb d = ζ 33 = ζ τ τ e = ζ ψ = ζ χ = 1 and all others equal to zero. Although we will take g Z = g U and m Z = m U / √ 2, these parameters do not have a significant influence on the dark matter abundance or LHC signatures considered (as long as g Z isn't too large and the Z isn't too light, as detailed below). Finally, in addition to the kinetic term for the leptoquark we also include non-minimal interactions between the leptoquark and the SM gauge fields appearing in gauge models, where T a = λ a /2 and λ a (a = 1, . . . , 8) are the Gell-Mann matrices.

Dark Matter Relic Surface
We now turn to the dark matter relic abundance. The presence of the coannihilating partner, ψ, can heavily impact the relic abundance of the dark matter candidate, χ 1 . We will restrict the parameter space by insisting that the model produces the observed relic abundance via thermal freeze-out. First, to clearly describe the physics of coannihilation, we will imagine that the masses of the dark sector (DS) particles χ 1 , χ 2 and ψ are above the electroweak scale and that the leptoquark and Z are heavy enough not to be present during χ 1 freeze-out (we allow them to be lighter in our numerical work below). If m ψ , m χ 2 m χ 1 then all processes connecting χ 1 to the thermal bath are heavily suppressed. As such, χ 1 will freeze-out when relativistic and will have a relic abundance many orders of magnitude too large to match observations (it will overclose the universe). However, the situation is very different if m χ 1 ∼ m ψ ∼ m χ 2 .
In this case the abundances of χ 1 , χ 2 and ψ are similarly Boltzmann suppressed during freeze-out of χ 1 . The conversion processes DS SM → DS SM have one rare, heavy particle and one bath particle in the initial state, so they have a rate exponentially larger than dark sector annihilation processes DS DS → SM SM, which requires two rare particles. In our model χ 1 χ 1 → SM SM is absent while both χ 1 ψ → SM SM (via an s-channel leptoquark) and ψψ → SM SM (via QCD processes) are efficient. This means that although χ 1 cannot efficiently annihilate with χ 1 , it can annihilate with ψ or be efficiently depleted by first converting into ψ, via χ 1 SM → ψ SM, which then annihilates via ψψ → SM SM. This coannihilation effect will be relevant if the mass splitting between ψ and χ 1 is small [59], and becomes more important as ∆ ψ shrinks. Similarly, if ∆ χ 2 0.3 then χ 1 χ 2 → SM SM via an s-channel Z may be efficient. However, since ∆ χ 2 2∆ ψ , g Z ∼ g U and m Z ∼ m U and since the relic abundance is predominantly determined by the most efficient process (which is χ 1 ψ → SM SM or ψψ → SM SM in the parameter space of interest), χ 2 does not play an important role in setting the relic abundance. We do however keep it in our numerical work.
To accurately calculate the relic abundance of χ 1 , the Boltzmann equation must track the abundances of χ 1 , χ 2 and ψ. This coupled differential equation can be written as a single differential equation [59], which is equivalent to the usual Boltzmann equation for a single species but with the annihilation cross-section replaced by where i, j ∈ {χ 1 , χ 2 , ψ} index the dark sector particles (note that ∆ χ 1 = 0), g i is the number of degrees of freedom of particle i (g χ 1 = g χ 2 = 2, g ψ = 4 in our model), the cross-section is σ ij = σ(ij → BB ) where B and B are particles in the thermal bath and x = m χ 1 /T . The effective number of degrees of freedom is given by We see that, for instance, when σ χ 1 ψ σ χ 1 χ 1 and both ∆ ψ and g eff are not too large, then Since Ω χ 1 h 2 ∼ 1/ σv , coannihilation effects then reduce the relic abundance of χ 1 .
To calculate the relic abundances we use model files written with FeynRules v2.3 [60] and solve the Boltzmann equations using micrOMEGAs v5 [61]. In principle, Sommerfeld corrections and the effects of bound state formation lead to corrections to the perturbative ���� ���� ���� ���� ���� Figure 1. The mass splitting ∆ ψ that yields the observed dark matter relic abundance as a function of the leptoquark mass m U and the dark matter mass m χ1 . While we set ∆ χ2 = 2∆ ψ , g Z = g U and m Z = m U / √ 2 in our numerics, the result is only weakly dependent on these parameters (as long as χ 1 χ 2 → SM SM via an s-channel Z is subdominant). Below the solid black line m χ1 < m U , while the dotted black line shows the resonant region where m χ1 + m ψ = m U . The observed dark matter relic abundance cannot be obtained in the white region.
cross-sections of coloured particle annihilation in the early universe. However, for the case of fermionic triplets in the TeV mass range, these corrections are negligible [62][63][64] and we do not include them. We use the relation between g U and m U given in eq. (2.20) and determine the value of the mass splitting ∆ ψ which will result in the observed relic abundance [65], The required value of ∆ ψ is shown in fig. 1. While we set ∆ χ 2 = 2∆ ψ , g Z = g U and m Z = m U / √ 2, we have checked that the result is only weakly dependent on these parameters. We see that the observed relic abundance can be obtained over a large mass range of χ 1 (100 GeV m χ 1 20 TeV) for mass splittings in the range 0.05 ∆ ψ 0.35. We do not consider leptoquark masses below 1.5 TeV due to collider constraints, or above 10 TeV, since g U becomes non-perturbative (due to eq. (2.20)).
We first discuss the region with m χ 1 < m U (below the solid black line). Along the resonant line m χ 1 + m ψ = m U the cross-section for the s-channel process χ 1 ψ → SM SM is large and so a large mass splitting is required to achieve the observed relic abundance. For a fixed m U , as m χ 1 reduces the process goes further off resonance, so a smaller mass splitting is required to compensate (to keep σ eff approximately constant). This continues until ψψ → SM SM becomes the dominant process at small m χ 1 . Once this happens the required mass splitting grows as m χ 1 , and so m ψ , reduces. In the region m U − m ψ < m χ 1 < m U , the mass splitting becomes smaller further from the resonance line, as the process again goes further off resonance.
When m χ 1 > m U (above the solid black line), in much of the parameter space the dominant process is χ 1 χ 1 → U 1 U 1 via a t-channel ψ (when m U m χ 1 , the leptoquark is still abundant as χ 1 freezes-out, so this process can efficiently annihilate χ 1 ). In this case the effective cross-section, eq. (3.2), is not exponentially sensitive to ∆ ψ , so small variations in m ψ cannot yield the observed relic abundance. This is indicated by the white region. Only when m U 3 TeV do the processes χ 1 ψ → BB and ψψ → BB , where B and B are bath particles, compete with χ 1 χ 1 → U 1 U 1 , so that ∆ ψ can be adjusted to obtain the observed relic abundance. Above m χ 1 ∼ 20 TeV the observed relic abundance cannot be obtained for any ∆ ψ or m U .

Coloured Coannihilation Partner Phenomenology
Since the dark matter candidate in our model is a Majorana fermion, the direct and indirect detection constraints are negligible [58,62,66]. As such, the strongest constraints on the simplified model parameter space will come from collider searches for the new particles. Collider searches for the leptoquark and Z have been well studied, and constrain the leptoquark to being heavier than ∼ 1.7 TeV [28,29]. Although the presence of the dark sector may weaken these limits, the leptoquark and Z will still have significant branching ratios to SM particles. The lightest dark sector particle, the dark matter candidate χ 1 , is a gauge singlet and only couples to the SM via the heavy leptoquark and Z , so searches for the dark matter candidate directly are challenging. However, as discussed in section 2, the coannihilation partner ψ is similar in mass to the dark matter candidate and is a colour triplet, making it directly accessible at hadron colliders. For this reason we focus on searches for ψ production.
After we fix the couplings as discussed in section 2.3, the model parameters relevant for collider searches are the masses of the dark matter candidate, m χ 1 , the coannihilation partner, m ψ , and the leptoquark, m U . In section 3 we found the ψ mass which leads to the observed dark matter relic abundance through thermal freeze-out, as shown in fig. 1. In our collider analysis we use this result to eliminate m ψ , allowing us to present our results in the m χ 1 -m U plane.
Since ψ is an unstable coloured fermion, it is important to compare its width with the QCD scale to see if it is likely to hadronize. We are primarily interested in the case where the Z 2 -odd fermions are significantly heavier than the SM particles, but lighter than the U 1 leptoquark (since the LHC will not be sensitive to dark sector particles heavier than a few TeV). Due to Z 2 -parity conservation, the dominant ψ decay channels are the threebody processes ψ → χ 1 bτ and ψ → χ 1 tν τ , which are mediated by an off-shell leptoquark. Due to the compressed spectrum, m ψ ∼ m χ 1 , the width of ψ is expected to be small. A straightforward computation (see eq. (B.3)) yields the partial width for the coannihilation regime ∆ ψ < 0.3 and for generic leptoquark couplings satisfying low energy constraints. 3 These suppressed widths for a TeV scale ψ suggest that it has enough time to hadronize into QCD bound states before it decays.

Ψ-Meson Spectroscopy
There are two types of mesonic bound states: 'psionium' states (ψψ) and 'open-psi' states (ψq), where q is a light SM quark. Their formation can be described, at leading order, by non-relativistic QCD using a modified coulomb potential similar to the hydrogen model [67], see appendix C for details. In this framework, the criteria for bound state formation is that the time it takes for the fermion pair to complete one rotation, t R , (defined in eq. (C.4)) must be larger than their intrinsic lifetimes, in our case τ ψ = Γ −1 ψ . In fig. 2 we show the regions of parameter space, on the relic surface, where the ψ lifetime is sufficiently long for the psionium states (ψψ) (denoted Ψ 0,1 , see below) and open-psi states (ψq) (denoted Ψ q ) to form. Note that the lifetimes of these bound states are not long enough for them to be considered stable on collider scales.
The coannihilating partner ψ thus gives rise to a rich spectrum of QCD bound states above the electroweak scale. The spectrum is expected to have a pattern similar to the bound state mesons containing charm or bottom quarks in the SM. These new heavy bound states can be classified using spin, parity, charge-conjugation and the new Z 2 -parity. We use the usual notation J P C ± with an additional ± subscript to indicate its Z 2 -parity. All psionia states are Z 2 -even, and the lightest psionium states are the S-wave pseudoscalar 0 −+ + followed by the vector 1 −− + (in analogy to the η c and J/ψ mesons in the charm sector of the SM, respectively). We name the ground state pseudoscalar state Ψ 0 and the lightest Type Production Decay

Collider Signatures
Given that ψ carries both colour and charge, it can be produced in large numbers at hadron colliders, mainly via QCD interactions. Single ψ production is forbidden by Z 2 -parity invariance, making pp → ψψ the main production mechanism at the LHC. At threshold, the pair of ψ particles will bind to form Ψ-mesons, leading to either single production of Ψ 0,1 states, or pair production of Ψ q states. In table 2 we show the main production and decay mechanisms for Ψ 0,1 and Ψ q at the LHC. Note that there is no QCD qq → Ψ 0,1 production or Ψ 0,1 → qq decay due to angular momentum conservation and since Ψ 1 is a colour singlet [68].

Search for Psionium in the Dilepton Channel
The cleanest channels for discovering the psionium states are the processes pp → Ψ 0 → γγ and pp → Ψ 1 → + − , which lead to a resonant peak in the diphoton and dilepton invariant mass spectra at the mass of the psionium, m Ψ 0,1 ≈ 2m ψ . Since we found that the dilepton decay channel of the vector meson pp → Ψ 1 → e + e − , µ + µ − (depicted in left panel of fig. 3) produces more stringent LHC limits, we focus on this process. 4 The partonic cross-section is given in the narrow-width approximation bŷ In appendix D we give the respective expressions for the partonic production cross-sections and the leading order decay rates. We used the package RunDec [70] to take higher order loop-corrections in the running of the strong coupling into account (which is important since the process is highly sensitive to the value of the strong coupling). Finally, we convolve the partonic cross-sections with the parton distribution function (PDF) set PDF4LHC15 nnlo mc [71][72][73]. The total cross-section is shown in fig. 4, as function of q q Z/γ  m ψ . Using the experimental upper bounds on the total cross-section for vector resonances given in ref. [74], we find that this dilepton search excludes m ψ < 280 GeV. Notice that this limit is independent of the exact values of the model parameters, as long as ψ hadronizes.

ψ Pair Production
We now turn to pair production of Ψ q bound states. These states predominantly decay through ψ → bτ χ 1 or ψ → tν τ χ 1 , since they are Z 2 -odd, leading to collider signatures bbτ τ χ 1 χ 1 , btτ ν τ χ 1 χ 1 or ttν τ ν τ χ 1 χ 1 . Due to the large top mass, ψ → χ 1 tν τ is typically kinematically forbidden for dark matter masses below 1 TeV, so we focus on ψ → χ 1 bτ decay channel (see the right panel of fig. 3). Since m ψ m q , the light quark acts as a spectator and we can assume that Ψ q and ψ share the same mass and decay width. Therefore, in a collider environment, we may neglect the presence of the light quark altogether. ATLAS and CMS have performed searches for pair-produced scalar leptoquarks decaying into the bbτ τ + E miss T final state at √ s = 13 TeV with an integrated luminosity of 36.1 fb −1 [53][54][55][56]. Unfortunately, these searches do not set strong limits on the model since they rely on large p T selection cuts for the reconstructed b-jets and τ -jets. The b-and τ -jet coming from ψ are expected to be relatively soft, due to the compressed spectrum m ψ m χ 1 , and therefore fail to pass the basic selection criteria for these searches. Indeed, when recasting [53][54][55][56] we found very small selection efficiencies, and these searches provided no constraints on our parameter space. In the next section we propose a dedicated search strategy for this type of compressed spectrum signatures at the LHC. However, the failure of the visible decay products of Ψ q to pass the typical LHC search cuts means that mono-X searches may be sensitive. The picture is that the soft (undetected) bτ pairs and hard missing transverse energy recoil against hard electroweak or QCD initial-state radiation (ISR) that does pass the selection criteria. In ref. [75], the process pp → QQ+j is investigated, where Q is a heavy coloured particle that decays into a dark matter candidate (and SM particles) and j is a hard (i.e., p j T > 250 GeV) ISR jet. The authors performed a recast of the ATLAS inclusive monojet search [76] and derive lower bounds on m Q for benchmarks of different spin and colour representations. The results obtained at NLO QCD for a fermionic top-partner can be readily applied to our scenario. It is worth noticing here that the study fixes the mass gap ∆m = m Q − m DM 40 GeV in order to avoid transitions to multi-jet + MET signatures. The relevant part of the parameter space in our setup where this condition can be satisfied is m ψ < 400 GeV, as can be seen in fig. 1. Coincidentally, the lower limit derived in [75] for the mass of the fermionic top-partner stands around the same value. Therefore, we conclude that inclusive monojet searches exclude ψ masses below 400 GeV at √ s = 13 TeV with L = 36.1 fb −1 of data.

Single ψ Production
It is possible to produce a single ψ in association with the dark matter particle χ 1 via pp → χ 1 ψτ (see table 2). Such production processes would directly probe the leptoquark coupling, which must be large to accommodate the B-physics anomalies. Unfortunately, these are 2 → 3 body transitions and are therefore heavily phase-space suppressed compared to the 2 → 2 processes. Moreover, the partonic initial state contains one b-quark, so this process is further suppressed by the small bottom PDF. While interesting, single ψ production modes are not expected to set competitive limits on the parameter space.

Dedicated Search for ψ Pairs
We propose a novel LHC search strategy for ψ pair production, specifically pp → ψψ → bbτ − τ + χ 1 χ 1 assuming a compressed spectrum scenario, m ψ m χ 1 . We focus on the experimental signature 2j b + τ h τ + E miss T , consisting of a pair of b-tagged jets (j b ), one τ -tagged jet coming from a hadronically decaying τ -lepton (τ h ), one light lepton e or µ coming from a leptonically decaying τ -lepton (τ ), and a large amount of missing transverse energy (E miss T ) mostly coming from the dark matter particles and τ -neutrinos. 5 Notice that this signature is similar to the SLT (single-lepton trigger) category described in [53]. The irreducible backgrounds are dominated by tt pair production decaying into t → bτ h ν and t → b ν, and sub-leading contributions from Drell-Yan (pp → τ + h τ − + jets) and diboson production (pp → ZZ → τ + h τ − j b j b ). The main reducible backgrounds consist of processes leading to objects misidentified as τ -jets. These fake taus arise from any process generating a prompt light lepton in association with jets and b-jets, with one of the QCD jets mistagged as a hadronic τ -lepton, j → τ h . In order to simulate the signal and background samples, we first implemented in FeynRules [60] the interaction Lagrangian eq. (2.18) and generated the UFO model [77] for the Monte Carlo event generator. The parton level event samples were simulated using MadGraph5 [78] and both showering and hadronization were performed in Pythia8 [79]. Jets were clustered using FastJet3 [80], while object reconstruction and detector effects were simulated with Delphes3 [81].
For the signal process, we simulated event samples for different points in the (m χ 1 , m U ) mass plane with the mass splitting ∆ ψ fixed by the observed relic abundance constraint, as shown in fig. 1. Moreover, the leptoquark couplings were fixed by eqs. (2.20) and (2.22) in order to satisfy the B-anomalies. For the background processes, in our analysis we only included the leading tt process given that the other irreducible backgrounds are subleading. Since the j → τ h mistagging efficiency are not easy to estimate without using data-driven techniques, we did not include these backgrounds in our analysis. Note that such backgrounds, while typically smaller than the leading tt process, are a non-negligible background that must be included in a more thorough analysis.
In the following we describe two analyses. One is a cut-based search using a single high-level observable R T (defined below) as a signal/background discriminator, while the other is a multi-variate analysis using a BDT classifier based on multiple low-level and high-level observables.

Event Selection and Key Observables
We now describe the basic event selections used in our analyses. Jets were clustered using the anti-k T algorithm with the standard narrow cone radius of R = 0.4 and candidate leptons = e, µ were selected if they passed the relative isolation requirements I rel = 0.05 within an isolation cone of radius R iso = 0.2. In order to improve the signal acceptance, it is important to relax the transverse momentum cuts as much as possible for any of the visible reconstructed objects since these are expected to be soft, due to the compressed spectrum. In our analysis we required p T (j) > 20 GeV, p T (τ h ) > 20 GeV, and p T ( ) > 5 GeV for jets, τ -jets and isolated leptons, respectively. The rapidity selections were taken to be the standard ones: |η(j)| < 2.5, |η(τ h )| < 2.3, and |η( )| < 2.5 with isolated electrons removed from the end-cap region (1.37 < |η(e)| < 1.52). For the b-tagging and τ -tagging efficiencies we used the default values of the ATLAS card in Delphes3. Finally, events are selected if 5 Signatures with two hadronic taus (2j b + 2τ h + E miss T ) can also be considered. While this channel has a larger branching ratio, it suffers from important reducible backgrounds (e.g., fake taus) that are difficult to determine using MC simulations. Estimating the LHC sensitivity of a search strategy for this signal category is therefore beyond the scope of this paper.  they contain two or more jets, of which at least one must be b-tagged, exactly one τ -jet and exactly one light lepton , such that the τ h pair has opposite electric charge.
After these basic selection cuts, we looked into different possible observables capable of distinguishing between the signal and the tt background. The two simplest discriminating quantities for this task are the missing transverse energy, E miss T , and the total amount of visible transverse momentum, defined as i.e., the scalar sum of the transverse momenta of the leading and sub-leading jets (with at least one being b-tagged), the τ -jet and the lepton. On one hand, because of the presence of the dark matter particles and neutrinos, signal events will tend to have more missing energy than tt events. On the other hand, because of the compressed spectrum, the visible energy in signal events will tend to be softer than in tt events. In the first two upper panels in fig. 5 we show these two normalised distribution for two signal benchmarks val-ues (m ψ , m χ 1 , m U ) = (1.12, 1.0, 1.5) TeV (red) and (m ψ , m χ 1 , m U ) = (0.45, 0.41, 1.5) TeV (blue) that satisfy the dark matter relic abundance requirement, and the tt background (shaded grey). In order to examine the relevance of angular information between final state objects, we also looked into the 2-point energy correlation functions. These observables are defined as for a set x 1 , · · · , x n of high-level reconstructed objects in the event. 6 Here d ij is the angular distance between two objects x i and x j , with transverse momenta p T i and p T j , respectively. We considered two possible distance functions: the azimuthal distance Φ ij ≡ 2(1 − cos ∆φ ij ) and the plane distance R ij ≡ (∆y ij ) 2 + (∆φ ij ) 2 , where ∆y ij = y i − y j , ∆φ ij = φ i − φ j are determined by the rapidity y i and azimuth angle φ i of the object x i . Notice that for d ij = Φ ij in eq. (5.2), each term reduces to the (squared) transverse mass M 2 T of a pair of objects, and the observable corresponds to the square of the total transverse mass M tot T of the group of objects x 1 , · · · , x n . If instead we fix d ij = R ij , then each term in eq. (5.2) corresponds to the k T -distance between the pairs of objects. We denote this observable k tot T and refer to it as the total k T -distance between x 1 , · · · , x n . In fig. 5 (left lower panel) we show as an example the k tot T observable computed for the set of visible final states τ h , ± , j 1 and j 2 , for signal and background. Notice that these observables containing pairwise angular correlations do not substantially improve the signal separation compared to observables such as S T without angular information (right upper panel). 7 M tot T , based on d ij = Φ ij , gives very similar results. Individually, none of the observables mentioned so far lead to a good separation between signal and background. Interestingly, much better discriminators are obtained by combining pairs of these observables. We found that the ratio between the invisible and visible total transverse momentum of the event, defined by considerably improves the signal separation compared to the individual E miss T and S T observables. This can be seen in the lower right panel of fig. 5, where we show the normalised R T distributions. Moreover, the shapes of the distribution for the signal process are independent of the particle masses m ψ , m χ 1 , which is not the case for E miss T . Finally, we also checked other ratios similar to eq. (5.3) by replacing the denominator with M tot T and k tot T . These ratios were found to have comparable performance to, and the same properties as, R T . We chose to cut on R T in our search, and we present the results in section 5.3.  Figure 6. Left panels: BDT scores for the two different benchmark points showing signal vs background discrimination. Right panels: ROC curves for the same benchmarks showing the tt background rejection rate vs signal efficiency for three different BDT classifiers based on only highlevel observables (purple curve), low-level observables (blue curve) and both high-level and low-level observables (green curve). For comparison, we also include the cut-based ROC curves for the R T ratio (dashed yellow) and E miss T (dashed red). Better performance is achieved by ROC curves stretching towards the upper-right corner and AUC values closer to 1.

Multivariate Analysis
Given that a combination of observables (R T ) performs better than the individual ones (E miss T and S T ), we also implemented a multivariate search using a BDT classifier. After imposing the same event selections as in section 5.1, a BDT classification score was extracted using low-level and high-level observables as input variables. The low-level observables consist of the p T and η distributions of the relevant reconstructed objects (j 1,2 , τ h and ) while the high-level observables consist of E miss T , S T , M tot T , k tot T and the R T ratio defined in eq. (5.3).
We tested three different BDT setups: (i) the high-level BDT trained with only highlevel observables, (ii) the low-level BDT trained with only low-level observables, and (iii) the combined BDT trained with both high-level and low-level observables. For each both signal benchmarks, we prepared an event sample with ∼ 45 K labelled events and a signal to background ratio of s/b ≈ 1.2. The samples were then split into 80% training and 20% testing sub-samples. The training was performed with XGBoost [83] using 100 trees with a maximum depth of 6, learning rate of η = 0.1 and a binary logistic for the objective function.
The resulting scores extracted from the combined BDT classifier, trained with the hyperparameters described above, is shown in the first column in fig. 5, for two signal benchmark mass points. A very clear separation between each benchmark signal and the tt background is obtained for this classifier. Furthermore, we measured the performance of the BDT classifiers using the receiving operating curves (ROC) displayed in the second column in fig. 6. There, one can observe that the three BDTs perform well (solid curves), achieving tt rejection rates above 10 for a signal efficiency of 50% and area-under-the-curve (AUC) values above ∼ 0.8. Overall, the low-level BDT (blue curve) performs better than the high-level BDT (purple curve), but the best classifier is the combined BDT (green curve) with AUCs above 0.9. The three BDTs considerably outperform the cut-based single discriminators (dashed curves). In this case the best cut-based observable was the R T ratio which produces moderate background rejection with an AUC close to 0.8.
These results demonstrate that a BDT classifier trained on both low-level and highlevel observables can improve the tt background rejection rate by at least a factor of two compared to the cut-based strategy based on the best high-level single observable, the R T ratio.

Results
A set of signal event samples for different values of (m χ 1 , m U ) were generated covering the parameter space region 200 GeV < m χ 1 < 1000 GeV and 1500 GeV < m U < 5000 GeV, with a rectangular grid with spacing (δm χ 1 , δm U ) = (50,200) GeV. In order to get statistically significant MC samples (especially in the tails of the R T distribution) we simulated around 1 M signal events for each mass point of the grid. For the background, we generated 500 K tt events.
For the cut-based strategy, after performing the event selections for each sample, we performed a statistical analysis using the binned R T distribution in the interval [0.5, 1.0] with a bin step of 0.1 including overflow in the last bin. Events were estimated with run-II integrated luminosity of L int = 140 fb −1 and a centre-of-mass energy of √ s = 13 TeV. In each R T bin, we fixed a 10% uncertainty for the tt background, which is approximately the same uncertainty that was estimated for the background S T distributions in ref. [54]. For the multivariate analysis, we trained the combined BDT on each point of the rectangular grid using the hyperparameters described above but with 50% training and 50% testing sub-samples (in order to retain enough MC events in the testing sample). We then used the binned score of the BDT classifier with 20 equidistant bins over the unit interval. Following ref. [54], we assumed a 10% background uncertainty for the BDT score. of the two search strategies, the expected exclusion limits were extracted at each point in the m χ 1 -m U plane using the asymptotic approximation of the CL s criteria [84,85] with the profile likelihood ratio as test statistic. This statistical model was implemented using the pyhf python package [86,87]. A mass point in the m χ 1 -m U plane is excluded at 95% confidence level (CL) if CL s < 0.05 is satisfied. In fig. 7 we present the final upper exclusion regions in the m χ 1 -m U plane given by the red area for the cut-based analysis and the purple area for the multivariate analysis. For comparison, we have also included the model-independent limits from the psionium production pp → Ψ 1 → + − (black dashed), as well as the expected reach of a monojet search (black dotted). We can see that both the cut-based and multivariate dedicated searches will give stronger limits that these recasted limits. Moreover, the multivariate analysis significantly outperforms the cut-based limits.

Conclusions
In this work we have investigated coannihilating dark matter, motivated by the vector leptoquark explanation of the B-physics anomalies. Assuming that the leptoquark is a gauge boson of a spontaneously broken gauge symmetry, and that dark matter is a fermion contained in a multiplet of that symmetry, the dark sector will also contain a coloured coannihilation partner. Furthermore, introducing UV-motivated dimension-5 operators leads to a Majorana dark matter candidate which is similar in mass to the coannihilation partner.
We determined the mass splitting between the Majorana dark matter particle χ 1 and its coannihilation partner ψ required to reproduce the observed dark matter relic abundance via thermal freeze-out (see fig. 1). Interestingly, processes with a leptoquark mediator are significantly more efficient than those with a Z mediator. We found that the relic abundance can be satisfied for dark matter masses up to ∼ 10 TeV and for leptoquark masses in the motivated range from 1.7-10 TeV. We then analysed the phenomenology of the model for the slice of parameter space which reproduces the observed relic abundance. Direct and indirect constraints are negligible as the dark matter candidate is a Majorana fermion, so collider experiments are the best probe of the parameter space.
Since leptoquark searches have been well studied and searches for the dark matter candidate in our setup are challenging, we focused on the coloured coannihilating partner, which can be pair produced with a relatively large cross-section at the LHC. The coannihilating partner lifetime is long enough for it to hadronize, forming both 'psionium' and 'open-psi' bound states. Psionium can decay via electroweak interactions into dileptons, leading to a lower bound of 280 GeV on the mass of the coannihilating partner. Open-psi predominantly decays to dark matter along with a b quark and a τ lepton via an off-shell U 1 leptoquark. Due to the compressed spectrum, the b quarks and τ leptons will be soft and may not pass selection cuts. We found that existing searches for bbτ τ + E miss T do not place constraints on the relevant parameter space. However, monojet searches are sensitive and currently exclude coloured coannihilation partner masses below 400 GeV.
We proposed a new search strategy centred around a new observable: the ratio between visible and missing transverse energy of the process. We performed both a cut-based and a multi-variate analysis that directly probe the mass of the coloured partner, assuming LHC run-II luminosities. Using the relic surface, we determined expected limits on the dark matter mass: m χ 1 400 GeV and m χ 1 600 GeV for the two analyses, respectively. We therefore conclude that a multi-variate analysis would significantly improve on current limits and could probe a significant portion of viable parameter space.
Our final results are summarised in fig. 7. Although this analysis was motivated by dark matter and the B-physics anomalies, we emphasise that (i) bound state formation and decay could be relevant for other coannihilating scenarios with coloured partners, and (ii) the ratio of visible to missing transverse energy could be a powerful discriminator in scenarios with compressed spectra and long-lived invisible particles.
Finally, while our results are encouraging, a more detailed analysis performed by the LHC experimental collaborations would provide a more robust estimation of the limits. For instance, including the 2j b + τ h τ h + E miss T signal category would strengthen the limits, and better accounting for mistagging rates and reducible and irreducible backgrounds would improve the accuracy of the limits. Furthermore, it would be necessary to carefully consider the trigger requirements for processes with soft final state objects, e.g., the possibility of triggering on our new observable R T could be investigated.

Acknowledgements
The authors would like to thank Ben Kilminster, for asking us questions that we hope this work goes some way to answering, Gino Isidori, for collaboration in the early stages

A 4321 models
We now briefly describe a SM extension that gives rise to a TeV scale U 1 ∼ (3, 1, 2/3) vector leptoquark as a gauge boson with non-universal couplings to quarks and leptons. These so-called 4321 models [20][21][22][23][24][25][26][27] are defined by the gauge group where SU (2) L is identified with SM weak isospin and both colour and hypercharge are diagonally embedded as SU The breaking G 4321 → G SM is typically induced by at least two scalar fields, Ω 1 ∼ (4, 1, 1, −1/2) and Ω 3 ∼ (4, 3, 1, 1/6), each developing a vev around the TeV scale. A Higgs doublet H ∼ (1, 1, 2, −1/2) is necessary for electroweak symmetry breaking. The resulting gauge sector (in addition to the SM gauge bosons) consists of a colour octet G ∼ (8, 1, 0), a colour singlet Z ∼ (1, 1, 0) and the leptoquark U 1 ∼ (3, 1, 2/3), all with masses above 1 TeV. There are two common 4321 implementations distinguished by the charge assignments of the SM matter fields: • (i) Standard 4321: All three would-be SM fermion families q i L , i L , u i R , d i R , e i R , are singlets under SU (4) and have the usual SM charges under SU (3) c ×SU (2) L ×U (1) Y . These fields therefore do not couple directly to any of the massive SU (4) gauge bosons. Three vector-like fermions Ψ i L,R ∼ (4, 1, 2, 0) are introduced in order to generate nonuniversal effective couplings between the SM left-handed fermions and the U 1 gauge leptoquark via fermion mixing, which arises from the Yukawa interactions¯ i L Ω 1 Ψ j R andq i L Ω 3 Ψ j R after spontaneous symmetry breaking. The matter content is shown in the first block for i = 1, 2, 3 in table 3.
• (ii) Flavoured 4321: In this case the would-be third family quarks and leptons are unified into SU (4) mulitplets Ψ 3 L = (q 3 L , 3 L ) T , Ψ 3+ R = (t R , ν R ) T and Ψ 3− R = (b R , τ R ) T and couple directly to the U 1 leptoquark. The first two would-be SM generations are SU (4) singlets and couplings to U 1 are induced via fermion mixing with two vectorlike fermions Ψ 1,2 L,R , as in the standard 4321. The matter content is shown in the first and second blocks for i = 1, 2 in table 3. Besides the minimal field content described above, more scalar or fermion fields are sometimes necessary to satisfy additional phenomenological requirements. For instance, symmetry breaking scalars Ω 15 transforming in the adjoint representation 15 of SU (4) can be included in order to induce mass splittings between gauge bosons and fermion components [24,25]. These fields will couple to SU (4) fundamentals asΨΩ 15 Ψ and a vev along the T 15 = diag(1, 1, 1, −3)/ √ 6 generator will lead to a mass splitting of order Ω 15 between the coloured and colourless components of Ψ. Other 4321 matter extensions require the presence of fermion singlets (1, 1, 1, 0) that give rise to light Majorana neutrinos through the inverse seesaw mechanism [22,88]. Finally, the 4321 models can be viewed as the low-energy limit of a more fundamental theory, such as the Pati-Salam cube model (PS 3 ) [23] (which can be embedded into a warped 5D construction [88,89]), the twin Pati-Salam model [27] or strongly coupled models with extended hypercolour [26].

B Partial Width Formulae
The tree-level partial decay width of the U 1 vector leptoquark into massive SM quarks and leptons is while the partial width into final state dark vector-like fermions is where we have used Källen's function λ(x, y, z) ≡ x 2 + y 2 + z 2 − 2xy − 2yz − 2zx.

C Bound State Formation
To describe the bound states we use non-relativistic QCD with a modified-hydrogenic model of a single-gluon exchange potential [67]. For coloured particles of mass m Λ QCD , the potential takes the Coulombic form, where C is a colour factor. For a colourless bound state C is simply the quadratic Casimir of the constituent particles (i.e., C = 4/3 for the psionium). We definē as the running strong coupling evaluated at the scale of the average distance r B between the two constituents. This distance is of the order of the Bohr radius a 0 = 1/(Cᾱ s µ), where µ is the reduced mass of the system. More precisely, for an S-wave ground state r B = √ 3a 0 . Bound states form if the revolution time, t R , is larger than the lifetime of the coloured particles. Combining the ground state energy with the virial theorem yields To determine the production cross-sections and decay widths of Ψ 1 we follow a generalisation of the results that apply to quarks and quarkonia [68]. The partial width for Ψ 1 decaying into a final state X can be written as Γ(Ψ 1 → X) → Γ(ψψ → X)|Ψ(0)| 2 , (D.1) where Ψ(0) is the wavefunction at origin for the ground state of the bound system and Γ(ψψ → X) is the corresponding process with a free pair of ψ andψ in the initial state.
In the modified-hydrogenic model, we obtain the following wavefunction where M is the psionium mass. The leading production mechanisms of Ψ 1 are: i) Electroweak production from qq. In the approximation that the psionium mass M is much larger than the Z boson mass, we may write the electroweak cross-section via an s-channel photon or Z boson as where x(x + 1) 2 + 4 ln x (x + 1) 3 .
(D. 6) Production in association with a photon or a Z boson is subleading and not considered here.
For m ψ ≈ 500 GeV and √ s = 13 TeV, the electroweak production constitutes almost 80% of the total cross-section, due to the running of the couplings and the phase-space suppression of gg → Ψ 1 g. To determine the cross-section eq. (4.2), we need to compute the branching fraction of the Ψ 1 bound state into the following states: i) SM fermions. The decay can proceed through a photon or a Z boson and the rate for fermions f L , f R is given by where n c = 1 for leptons and 3 for quarks.
ii) ggg or γgg. The decay rates to three gauge bosons (which leads to the bulk of the hadronic decay modes) are Γ Ψ 1 →ggg = 5 π 2 − 9 27π Assuming that there are no other decays with considerable rates, the branching ratio to any single flavor of leptons is around 10%.