Non-perturbative effects for dark sectors with QCD portals

In this work, we consider a class of dark matter (DM) models where the DM does not directly interact with the Standard Model (SM) particles at the tree-level. Therefore, the coannihilation mechanism is crucial in achieving the correct DM relic abundance, which in turn requires the coannihilating partner to be close in mass to the actual DM particle. In our systematisation of the models' class, the mediator and the coannihilation partner are assumed to be charged under QCD interactions. This last feature calls for a scrutiny of non-perturbative effects, namely Sommerfeld factors and bound-state formation, on the annihilations of the colored partner. Such non-perturbative effects are illustrated with an example model comprising a scalar leptoquark mediator, a Dirac vector-like fermion coannihilation partner, and a singlet DM fermion. Phenomenological features of this model, namely DM direct and indirect detection prospects, collider implications, and impact on the muon anomalous magnetic moment, are discussed.


Introduction
The dark matter comprises about eighty-five percent of the matter in the universe.However, the Standard Model (SM) -the most successful theory in particle physics -fails to provide a dark matter candidate.Therefore, going beyond the SM (BSM) appears to be inevitable.In the context of particle physics, dark matter can be made of one or more new particles (see e.g.[1,2] for extensive reviews).Dark matter particles must be stable on cosmological time scales and expected to be uncolored, electrically neutral, and weakly interacting.Symmetries beyond the SM are typically required to stabilize a dark matter candidate.In the standard thermal freeze-out scenario, dark matter particles were in thermal equilibrium in the early universe, and later on they annihilated into particles of the visible sector.Today, in the universe, we observe the relics, namely the leftover dark matter abundance that has survived the annihilations as of now.The relic dark matter energy density is a precisely measured cosmological quantity, and the Planck collaboration provides Ω DM h 2 = 0.120 ± 0.001 [3].
There have been a plethora of particle physics models explaining the origin of dark matter.Depending on the model details, the dark matter may or may not couple directly to the visible sector.Due to increasingly stringent experimental constraints [4], there has been renovated interest in dark matter candidates that are very weakly coupled to the SM sector.This can be realized in various ways, which include feebly interacting massive particles (FIMP) [5][6][7] or gravitational dark matter [8][9][10][11][12][13][14] (in the latter case the only interaction between the hidden and visible sector is mediated by gravity).
In this work, we consider a class of models where the dark matter does not have any direct interaction with the SM sector.As a result, the dark matter pair annihilation crosssection to the visible sector is typically small.In order to achieve the correct dark matter relic abundance, an additional dark partner with a large coannihilation cross-section is then often required.Such a large cross section is easily attainable if the coannihilation partner carries SM charges, and then has sizable couplings to the SM particles.For a scenario of this type, three different sectors are needed: (i) the visible sector, (ii) the dark sector, and (ii) the mediator sector.In our framework, the dark sector consists of a Majorana or a Dirac fermion dark matter, χ, which is a singlet under the SM.In addition, the dark sector contains a coannihilation partner, which is a dark Dirac fermion, ψ, that transforms non-trivially under the SM.As for the mediator sector, that couples to both the SM and dark sectors, we assume it is made of a scalar particle, φ.In particular, we focus on colored mediators and coannihilation partners with non-zero hypercharge, even though they can also be charged under the weak isospin.
In this setup, the dark matter elastic scattering with the SM particles take place only via loops, hence the non-observation of dark matter signals in the direct detection experiments can be naturally explained while retaining a rich phenomenology at collider facilities.The presence of a colored scalar φ, in particular a leptoquark option, may also have an interplay with the anomalous magnetic moment of the muon.To stabilize the dark matter, we impose a Z N symmetry under which all dark sector particles, namely, χ and ψ are charged.
Since the mediator, as well as the coannihilation partner, carry color charges, a standard estimation of the dark matter relic abundance in the coannihilation scenario is not accurate.More specifically, a precise calculation involving non-perturbative effects, namely the Sommerfeld factors [15][16][17] and bound-state formation [18,19], must be considered for the pair annihilations of the colored coannihilation partners (we refer to them as nonperturbative or near-threshold effects along the paper).The two effects are the manifestation of multiple soft exchanges of QCD gluons and, therefore, they should be both included in the relevant cross sections.Recent investigations [16][17][18][20][21][22][23][24][25][26][27][28][29][30][31][32][33][34] have shown that the mass benchmarks that give DM energy densities compatible with cosmological observations are rather different from the case with no threshold effects.Accordingly one finds important changes of the model parameter space that is compatible with the observed DM energy density.
In this work, our particular focus is to precisely extract the dark matter relic abundance by accounting for the nature of the colored coannihilation partners and mediators of our framework.To this aim, we consider a particular exemplary model where the DM is accompanied by a color-carrying coannihilation partner and a scalar leptoquark, which has been formerly considered in the literature [35] (see also [36][37][38]).The importance of coannihilations of the colored coannihilation partner has been highlighted in order to obtain the correct DM energy density.With respect to earlier studies, we include the Sommerfeld factors, as well as the bound-state formation and decay, for the coannihilating partner, namely, the dark vector-like fermion (a Dirac fermion that carries a dark charge and is vector-like under the SM gauge group).The nature of the dark matter fermion dictates the relevant pair annihilations of the dark vector-like fermion.For the Majorana DM option, particle-particle (ψψ) annihilations are possible in addition to particle-antiparticle (ψ ψ) annihilations.Vector-like particle-antiparticle pairs combine either in color-singlet or color-octet states, whereas ψψ pairs organize in color antitriplets and color sextets.Color-singlet and antitriplet configurations feature an attractive potential and can sustain bound states.We compute the corresponding bound-state formation cross section and bound-state decays in the framework of non-relativistic effective field theories (NREFTs), and include their effect in the numerical extraction of the DM energy density.
The paper is organized as follows.In Sec. 2, we introduce the framework for our class of models.In Sec. 3, the dark matter relic abundance is computed by taking into account non-perturbative effects.Dark matter direct and indirect detection prospects, and collider constraints on the model parameters, as well as correlated phenomenologies, are summarized in Sec. 4. Finally, conclusions and outlook are given in Sec. 5, whereas supplementing material is collected in the appendices.

Model Setup
In this section, we discuss the construction of the models' class and specify the relevant interactions of the dark sector with the visible sector.We first discuss the general framework and then focus on an exemplary model.

General framework
In this work, we consider a scenario where the dark matter, χ, is a gauge singlet under the SM group.The dark matter talks to the SM sector via a colored scalar mediator, φ, and a colored fermion, ψ, which is the coannihilation partner.To cancel the gauge anomalies, we consider the fermion ψ to be vector-like under the SM and refer to it as a dark vector-like fermion (DVLF), since it belongs to the dark sector.To form a gauge as well as Lorentz invariant interaction of the mediator with the dark sector, the BSM scalar and the fermion must carry the same quantum numbers under the SM group.
particle type dark charge interaction with SM dark interaction    Dark matter annihilation χχ → SM SM (χχ → SM SM SM SM) takes place only (already) at the loop-level (tree-level).For the coannihilation partner, both the ψψ → SM SM and ψψ → SM SM SM SM occur at the tree-level.Moreover, coannihilation channels χψ into visible sector particles lead to SM SM and SM SM SM final states at the tree-level.Processes relevant for dark matter direct detection (for example, from bottom to top in the leftmost diagram for χ) only happen at loop-level.
In table 1, F and S stand for fermion and scalar respectively.Moreover, a SM gauge boson is denoted by V. Note that if the scalar φ does not couple to SM fermion bilinears, then χχ annihilation to SM sector occurs only at the loop-level.For example, χχ → gg can occur at one-loop level, where φ and ψ states propagate inside the loop4 .However, in a general scenario, dark matter pair annihilation to the visible sector is allowed at the treelevel.Moreover, pair annihilation of the DVLF as well as dark matter-DVLF coannihilation processes already take place at the tree level.Annihilation and coannihilation channels of the dark sector particles are schematically presented in figure 1.

An example model
In this work, we shall explore a specific model realization of the more general framework presented in section 2.1.In order to assess the relevance of non-perturbative effects, and highlight possible experimental signature and constraints, we consider φ to be a LQ that carries hypercharge Y φ = 1/3 and is a singlet under SU(2) L .Although we scrutinize a particular model realization, we acknowledge that each model merits an independent study on its own.Moreover, we study two scenarios where the dark matter is a (i) Majorana fermion and (ii) Dirac fermion.For the former case, in order to achieve the stability of the dark matter, we assign a Z 2 odd charge to only χ and ψ, whereas the rest of the particles are even.The full quantum numbers of the BSM fields under the SM×Z 2 are as follows: On the other hand, for a Dirac fermion dark matter scenario, a Z 3 symmetry is imposed with ω 3 = 1.The dark matter interaction to the DVLF and the mediator, which induces DM annihilations to the visible sector, is realized by the following interaction: where we have used the usual notation χ c = χ T C, where C is the particle-antiparticle conjugation matrix.And for the Dirac case, we have, The coupling y can be made real in both the cases by a field redefinition.
The LQ is neutral under Z N and has direct interactions with the SM fermions [39], Here, i, j = e, µ, τ are the family indices.We assume baryon number conservation by assigning the LQ and the DVLF the baryonic charge B = −1/3 under the global U (1) B .Consequently, diquark couplings of the LQ are absent, and the theory is safe from proton decay.The theory may also have a global lepton number symmetry U (1) L as in the SM.Therefore, φ can carry a negative unit of lepton number; consequently, we assign the same lepton number to ψ.This, however, does not lead to any phenomenological implications 5 .The LQ couplings y L ij and y R ij to the SM fermions are a priori free parameters.However, these couplings are constrained from the direct searches at the colliders as well as from flavor violating processes, which will be discussed later in the paper.In the next section 3, when discussing the dark matter relic abundance, we turn on only a single coupling from each y L,R and take them to be equal to one.At this point, specifying which non-zero entry is taken, however, is irrelevant.In Section 4, where we present phenomenological implications of the model, we explicitly specify the texture of the relevant Yukawa couplings.

Dark matter relic density
The dark matter cosmological abundance is accurately determined by measuring the CMB anisotropies and it amounts to Ω DM h 2 = 0.120 ± 0.001 [3], where h is the reduced Hubble constant.It stands as the main observable that any compelling dark matter model has to comply with.Upon selecting a viable mechanism to produce dark matter particles in the early universe, one can use the observed relic density as a powerful constraint on the model parameters.
In this work we consider thermal freeze-out [40,41].For such a mechanism to work, dark matter particles have to be kept in equilibrium at high temperatures, and the relevant processes for determining the relic density are dark-matter pair annihilations.When the temperature of the expanding universe drops below the dark matter mass, the corresponding particle densities become Boltzmann suppressed and the annihilations cannot keep up with the expansion of the universe.The chemical freeze-out occurs at temperatures T ≈ m χ /25, therefore, dark matter particles are non-relativistic.
The presence of additional dark-sector particles during freeze-out may severely affect the relic density: one has to track the (co)annihilations of additional partners when these are close in mass with the actual dark matter particle [41,42].For the model under consideration, the DVLF plays the role of the coannihilating state. 6The thermal freezeout has been studied for this model in refs.[35,36], where it has been highlighted the role of coannihilations in order not to overclose the universe and relative mass splittings as small as ∆m ≡ (m ψ − m χ )/m χ ≈ 10 −3 are considered (see also [37] for the extreme case ∆m = 0).The impact of coannihilating processes depend strongly on (i) the mass splitting between the DM particle (χ) and the coannihilating species (ψ); (ii) conversion rates between the dark matter and the coannihilating partner that put them in thermal contact.For the portal coupling y, we consider a range y ∈ [0.01, 2], which well ensures fast conversion rates [43,44].A complementary dark-matter production mechanism for the model under study, namely the conversion-driven freeze-out (or co-scattering) [43,44], has been addressed in ref. [36], where much smaller Yukawa couplings 10 −8 y 10 −4 are considered.For even smaller portal couplings, the freeze-in mechanism is the viable option [5,6], and we leave it for future work on the subject.
In summary, in the coannihilation regime, non-relativistic DVLFs pair annihilation also contribute to the depletion of the dark matter.At variance with the dark matter particle, DVLFs feel QCD strong as well as electroweak interactions.Pair annihilations of slowly moving charged particles get affected by long-range interactions mediated by soft gaugeboson exchange, that induce near-threshold effects, most notably Sommerfeld and boundstate effects [15,16,18,19].Due to the observed hierarchy between the corresponding SM gauge couplings, strong interactions are largely dominant, and we focus on them in this work. 7

Boltzmann equation and cross sections
The effect of a co-annihilating partner (ψ) in thermal equilibrium with the actual dark matter particle (χ) can be captured by a single Boltzmann equation [40][41][42] where H is the Hubble rate of the expanding universe and n denotes the total number density of dark sector states χ and ψ.By assuming the dark matter being a Majorana fermion, the total equilibrium number density, which accounts for both the particle and antiparticle species of the dark sector is where p ≡ d 3 p/(2π) 3 , E p = m χ + p 2 /(2m χ ), p ≡ |p|, g χ = 2 and g ψ = 2N are the particles degree of freedom, i.e. spin polarizations and color N = 3, and the effective thermally averaged annihilation cross section reads [42] χχ pair annihilation:-The cross section of the contributing processes have to be handled with care for several reasons.A first observation is about the relative importance of χχ pair annihilations into SM particles.This process can occur via loops into a two-body final state, χχ → 2 SM, or via off-shell decays of leptoquarks into a four-particle final state χχ → φ * φ † * → 4 SM.As noted in [35], the first class of processes feature typical loopsuppression factors, whereas the latter are phase-space suppressed.Only for ∆m/m χ 0.4 [35], the loop-induced and off-shell decays of the leptoquark is comparable with the 2 → 2 coannihilation channels χψ → 2 SM and ψ ψ → 2 SM.In this work we restrict to smaller mass splittings, namely ∆m/m χ 0.1, as motivated by the need efficient coannihilation regime [35,36], and we find that the two classess χχ → 2 SM and χχ → φ * φ † * → 4 SM are negligible with respect to χψ → 2 SM and ψ ψ → 2 SM (in agreement with the estimations of ref. [35]).
The phase-space suppression of dark matter annihilations into four-body final state is lifted for m χ > m φ .In this case, the leptoquark can be produced on shell and, by relying on the narrow-width approximation (NWA) [45][46][47], the cross section may be factorized for a given final state in the form σ(χχ (here for the decay of φ and φ † in the same final state).Upon looking at the inclusive annihilation process, namely the sum of all possible four-body SM final states, the relevant cross section reduces to σ(χχ → φφ † ). 8 In order to properly use the NWA, there are several conditions to be fulfilled [47].As far as we are concerned, Γ φ /m φ 10 −2 1 for the y L , y R values adopted in this work (cfr.eq.(3.4)); the daughter particles are much lighter than the leptoquark; the leptoquark propagator is separable.However, at the opening of the channel, namely m χ m φ , we are not far from the mass threshold, which is another condition for using safely the NWA to be fulfilled. 9Moreover, since annihilations happen in a thermal environment, the so-called forbidden region opens up for m χ slightly smaller than m φ , the negative mass gap being compensated by the thermal kinetic energy of the incoming dark matter pair [41].We take into account the forbidden region when computing the thermally averaged cross section σ χχ v rel .Both in the forbidden region and the allowed region with (m χ − m φ )/m χ 0.1, the velocity expansion cannot be performed [41], and we use the exact expression for the cross section of the processes χχ → φφ † and χ χ → φφ † (Majorana and Dirac DM option respectively).Away from the threshold region, we find that σ χχ v rel has a leading p-wave contribution for the Majorana case, whereas Dirac dark matter annihilation σ χ χv rel feature an s-wave leading contribution.The exact cross sections for the processes χχ → φφ † and χ χ → φφ † , as well as the velocity-expanded ones, are provided in appendix A.
χψ coannihilation:-Next, the coannihilation process χψ → q c ¯ and χψ → Q c L proceed via an s-channel exchange of a leptoquark.The cross section become resonantly enhanced whenever m χ +m ψ ≈ m φ , which in the coannihilation scenario with small relative mass splittings gives m χ , m ψ ≈ m φ /2.The total decay width of the leptoquark has to be included to properly regulate the annihilation process. 10For order-one Yukawa couplings y, y L and y R , χψ → q c ¯ and χψ → Q c L give large cross sections in the region m χ ≈ m φ /2, that in turn produce prominent effects in the DM energy density (see figures 7 and 9).There is another coannihilation process, namely χψ → φ g → 3 SM, that proceed both via a t-channel exchange of the DVLF as well as an s-channel mediated by the leptoquark.The three-body final state (a gluon, a quark and a lepton) is obtained after the decay of the unstable leptoquark.Analogous arguments as discussed for χχ pair annihilation on the applicability/approximations of the NWA hold.The leptoquark decay width at leading order reads The first two contributions stem for the leptoquark decays into right-handed and lefthanded quark and lepton pairs respectively (we treat the SM particles as massless since in this work we take m φ of order 1 TeV).The third contribution appears only if m φ > m χ +m ψ .ψψ and ψψ pair annihilations:-Finally, the pair annihilation of non-relativistic DVLFs, both σ ψ ψv rel and σ ψψ v rel , can be affected by non-perturbative effects due to repeated gluon exchange, see representative diagrams in figure 2. The latter cross section is triggered only by Yukawa interactions and for the Majorana dark matter option via the process ψψ → φφ; see rightmost diagram in figure 2. For incoming scattering states, longrange interactions induce Sommerfeld factors, which enhance (suppress) the annihilations for an attractive (repulsive) potential experienced by DVLFs pairs.Moreover, there is an additional manifestation of repeated soft gauge-boson exchange: the presence of metastable bound states.The bound-state formation process, and the subsequent bound-state decay, triggers an efficient way to deplete further the QCD-charged coannihilating states, and then the overall dark matter abundance [29][30][31][32]50].The inclusion of bound-states can be implemented in the annihilation cross section of DVLFs in (3.3) through the following effective cross section [29,31] The corresponding effective cross sections for DVLF particle-particle and antiparticleantiparticle pair annihilation have the same form as eq.(3.5) with the subscripts in the cross sections and widths replaced by ψψ and ψ ψ.In eq.(3.5) the first term stems for the annihilation of scattering (or unbound) states, whereas the second term encodes the reprocessing of an unbound pair into a bound state.Here, the quantities that enter are the thermally averaged bound-state formation cross section σ n ψ ψ,bsf v rel , the decay width of the bound states Γ n ψ ψ, and the bound-state dissociation Γ n ψ ψ,bsd .The combination of the bound-state decay and dissociation widths Γ n ψ ψ/(Γ n ψ ψ + Γ n ψ ψ,bsd ) takes into account the ionization of a given bound state in the thermal environment, and dictates how efficiently the bound-state formation contribute to the depletion of colored pairs.Upon the inclusions of bound-to-bound transitions, eq.(3.5) has to be modified [30,51] and bound-state effects become even more relevant.As for the DVLF annihilations into leptoquark pairs, which mediate the subsequent decays into SM four-body states, we find that ψ ψ → 4 SM and ψψ → 4 SM are phase-space suppressed and practically negligible for m ψ < m φ with respect to ψ ψ → 2 SM (for the DVLFs the 2 → 2 annihilations occur via QCD interactions without the need of the φ mediator, see e.g.leftmost diagram in figure 2) .At the opening of the on-shell region for the leptoquarks, and in full analogy with χχ annihilations, the suppression is lifted and we approximate the inclusive cross section with We have checked our analytical expressions for the cross sections of 2 → 2 annihilation processes σ ij that enter eq.(3.3) with the model implementation in MadGraph [48], and the corresponding relic density with micrOmegas [52].However, in order to go beyond the free annihilations and include non-perturbative effects, an estimation of Sommerfeld and bound-state effects for colored DVLFs is needed.This is the subject of the following sections.

Near-threshold effects in NREFTs
The dark matter is a SM gauge singlet, hence the free cross section accurately account for the corresponding pair annihilation.Conversely, annihilating DVLF pairs, either ψ ψ, ψψ and ψ ψ annihilations, are affected by soft gluon exchanges (see figure 2 for exemplary diagrams).In the following, we assemble existing results, and obtain new ones, in order to compute relevant cross sections and decay widths of non-relativistic DVLF pairs in the early universe thermal environment.We exploit the hierarchy of energy scales that is typical for non-relativistic particles moving with relative velocity v rel , namely m ψ m ψ v rel m ψ v 2 rel , by replacing the fundamental DM theory with a sequence of non-relativistic effective field theories (NREFTs).For Coulombic bound states the relative velocity is fixed by the virial theorem as v rel ∼ α s , hence the corresponding hierchy of scales is m ψ m ψ α s m ψ α 2 s .In particular, we shall exploit the framework of NRQCD [53,54] and potential NRQCD (pNRQCD) [55,56], since the DVLF well qualifies as a heavy quark from the QCD perspective.The original formulation of such EFTs were conceived for heavy quark and antiquark pairs, then leading to color-singlet and color octet states.In our work, it will be relevant to address DVLF particle-particle pairs as well.We shall rely on the corresponding NRQCD and pNRQCD for two heavy colored particles (or antiparticles) as detailed in ref. [57].
pNREFTs are useful for our scope since they stand for the quantum-field theories of non-relativistic interacting pairs, both for scattering or bound states, and allow to systematically describe pair annihilations and pair-to-pair transitions.Since the relevant processes occur in the early universe, we exploit the formulation of pNRQCD at finite temperature [58][59][60].There has been a recent effort in transferring and adapting the NREFTs for dark matter freeze-out [24,25,30,32,33,50,[61][62][63].A detailed inspection of the interplay with thermal scales in the construction of NREFTs relevant for dark matter annihilations has been recently carried out in ref. [64].In this work, we restrict to the bound-state formation process as induced by the radiative emission of a gluon (see Sec. 3.2.2).Moreover, we discuss the applicability of NREFTs when the dark sector particles annihilate into final states with comparable masses (in the model at hand this means DVLF annihilations into leptoquark pairs).The latter situation invalidates the velocity expansion, and hence, some care is needed when one aims to include non-perturbative effects close to the opening of a mass threshold.

Sommerfeld factors for pair annihilations
The annihilation process of fermion-antifermion pairs is encoded in the imaginary part of the matching coefficients of four-fermion operators, that are organised according to spin and color representations [54].In pNRQCD, this translates into an imaginary local potential for the pairs, which is inherited from four-fermion operators of NRQCD [55,56,65].Whenever we consider DVLFs directly annihilating into light Standard Model particles, namely quarks and gluons, the large mass gap between initial and final states makes the wavelength of final-state particles of order 1/m ψ (and the energy scale for such annihilation being of order m ψ ).This scale is much smaller than the corresponding wavelength of incoming non-relativistic DVLF states, i.e. 1/m ψ v rel .For such a reason heavy-pair annihilations are well described by local interactions, i.e. the four-fermion effective operators of NRQCD (see figure 3), and the corresponding imaginary local potential in pNRQCD.
The factorization of hard modes and soft scales is a built-in feature of NREFTs [54][55][56].Soft gluon exchanges, which correspond to energy modes of order m ψ α s , are encoded in the real part of the potentials of color-singlet and color-octet pairs, which at leading order read11 where C F = (N 2 −1)/(2N ).The fermion-antifermion wavefunctions in pNRQCD, which at leading order in the multipole expansion are the solution of the Schrödinger equation with the potentials in eq.(3.6), accounts by construction for the effect of multiple soft gluon rescattering [55,56].Then, by combining the known results for the matching coefficients of heavy quark-antiquark annihilations [54], we can obtain the annihilation cross section for ψ ψ pairs into Standard Model QCD states, namely gluons and quarks, the latter counted by the number of flavors n f .The main advantage over exploiting the NRQCD framework is a transparent organization of the contributing partial waves, color and spin states to the DVLFs annihilation, which makes manifest the corresponding Sommerfeld factors.In the following, we provide the analytical expressions of the Sommerfeld corrected cross sections at leading order in the velocity expansion.This corresponds to the inclusion of the leading dimension-six operators of NRQCD.
The color and spin-averaged cross section that accounts for ψ ψ annihilation into SM gluons and quarks, (σ where the strong coupling constant that appears in the NRQCD matching coefficients is evaluated at the hard annihilation scale p (0)| 2 correspond to the squared wave function of the colorsinglet and octet pairs evaluated at the origin, because of the local nature of the annihilation process into light states.The Sommerfeld factors encode the soft contribution to pair annihilations.The color-singlet and color-octet Sommerfeld factors, with orbital angular momentum = 0, read ) where ζ = α s (µ s )/v rel .Here the strong coupling constant is evaluated at soft scale, µ s ≡ m ψ α s , namely the energy/momentum scale typical of soft-gluon exchanges. 12In eq.(3.7), DVLF pairs in a color-singlet annihilate into gluons only, whereas color-octet pairs can annihilate into both gluons and quarks, as one may see from the appearance of n f .
Besides direct annihilations into light two-body SM states, there are two additional processes which are driven by the Yukawa-portal and QCD interactions, namely ψ ψ → φφ † and ψψ → φφ and ψ ψ → φ † φ † , that mediate the DVLFs annihilation into a four-body final state via on-shell decays of the leptoquark pairs.For these processes, the velocity expansion breaks down when the DVLF and leptoquark masses are nearly degenerate (this has been highlighted in the context of dark matter freeze-out [41]).This is because, for small mass splittings ∆m ψφ ≡ m ψ − m φ , higher powers of non-relativistic velocity of the incoming DVLFs become of comparable size with ∆m ψφ , and it is not sufficient to retain the leading term in the velocity expansion of the cross section. 13In order to recast this situation in the language of effective field theories, let us take the corresponding annihilation diagram into leptoquark pairs, see figure 3 (right, diagram b).The typical momentum of the massive leptoquarks is parametrically of order |k φ | ≈ m ψ v rel , m ψ ∆m ψφ /m ψ , which qualifies as a small energy scale with respect to m ψ .Hence, the annihilations of DVLF into leptoquarks with slightly smaller masses cannot be described by local annihilations, or equivalently, by local four-fermion operators of NRQCD.The incoming DVLF particles do not have to come very close to annihilate, since the wavelength of the final-state particles is comparable 12 We evolve the strong coupling constant at one loop with the additional colored states φ and ψ, a scalar and a fermion respectively, as follows ∂tg2 s = , where t = ln μ2 parameterizes the MS renormalisation scale, NG = 3 indicates the SM generations, NF = 1 and NS = 1 stand for the BSM colored fields, and N = 3 for the specific case of QCD. 13 In the center of mass of the collisions, the final state particle momentum can be written as follows Upon expressing the denominator with the relative velocity as a geometric series and using the definition of the mass splitting ∆m ψφ = m ψ − m φ , one obtains for the final-state momentum with that of the incoming non-relativistic states.As a result, Sommerfeld and bound-state effects on pair annihilations in this regime, are expected to be less relevant.Lacking of a quantitative assessment for this situation, we do not include any of them at the opening of the mass threshold. 14  In practice, we assess the convergence of velocity-expanded cross sections to the exact cross sections as a function of the mass ratio m φ /m ψ .We have checked that for m φ /m ψ ≤ 0.8 the velocity expansion can be used, 15 and we include the Sommerfeld effect accordingly.This is also the range that makes the NWA performing rather good because we are sufficiently away from the leptoquark pair threshold, hence, the cross section of the process ψ ψ → φφ † describes well the annihilations into all possible four-body SM states.At leading order in the velocity expansion, we find that ψ ψ → φφ † contribute to the dimension-six spin-triplet color-singlet and spin-triplet color-octet operators of NRQCD.The corresponding spin-and color-averaged cross section reads (3.9) For Majorana fermion dark matter there is additional annihilation channel for DVLF pairs, namely ψψ → φφ and the complex conjugate process.In our assignation of the quantum numbers, ψ has the same SU(3)-color charge of a SM quark, cfr eq.(2.1).Hence, annihilating DVLF antiparticle-antiparticle pairs organize either in a color antitriplet or color sextet state, 3 ⊗ 3 = 3 ⊕ 6, whereas particle-particle pairs into the corresponding conjugate representations.In the following we simply denote the representations of the ψψ and ψ ψ pairs with [3] and [6], since a color antitriplet (sextet) has the same symmetry property of a color triplet (antisextet).The corresponding NRQCD and pNRQCD for identical fermions can be found in ref. [57].A color triplet pair feels an attractive potential, whereas the color sextet a repulsive one.The leading order potentials read 14 A qualitative estimate of the Sommerfeld and bound-state effects can be inferred by taking the Coulombic wave functions evaluated at a typical soft scale of order m ψ αs or m ψ v rel , rather than at the origin.S-wave bound-state decay widths would be then suppressed by a factor ∼ 1/e 2 .As for the scattering state wave function, the modulus squared of the corresponding hypergeometric function would enter.We checked that, for the attractive channels |Ψ  p (r)| 2 , a reduction of an order of magnitude is found with respect to their value at the origin when one insert r = 1/m ψ αs instead. 15Our finding compares well with the original statements in ref. [41], which considered the velocity expansion to be valid for mass rations of 0.85-0.9.Here we adopt a slightly more conservative condition. [1] [3] [8] [6] The spin-and colored averaged cross section is16 and the corresponding attractive Coulombic Sommerfeld factor reads with C a ≡ (N + 1)/2N .We notice that, at variance with the annihilation processes ψ ψ → gg, ψ ψ → q q and ψ ψ → φφ † , only the attractive color-antitriplet channel contribute at leading order in the velocity expansion.In agreement with general arguments on the symmetry of identical particle annihilations as given e.g. in ref. [66], we find that the velocity-independent cross section (3.11) corresponds to DVLF pair in a spin triplet.
In figure 4, we show the thermally average Sommerfeld factors for the different color representations.The thermal average is performed in the standard way, see e.g.[19], that amounts at a taking Maxwell Botlzmann distribution for the incoming DVLFs.The antitriplet and singlet Sommerfeld factors enhance the corresponding contributions in the free cross section, whereas the sextet and octet suppress them.The mass splitting is fixed to ∆m/m χ = 10 −2 , though we find that there is no appreciable difference for ∆m/m χ ∈ [10 −3 , 10 −1 ].

Bound-state formation, dissociation and decays
Bound-state formation is yet another manifestation of repeated soft-gluon exchanges: in the spectrum of a two-particle system there is an above-threshold continuum of states along with bound states below threshold.In this section, we address the bound-state formation, bound-state decays and bound-state dissociation processes.The latter is a genuine thermal process that happens in a thermal environment and describe the thermal break up of a bound state when hit by a sufficiently energetic thermal gluon.Its interplay with the bound-state decay dictates how efficiently DVLFs are depleted in the form of bound states.All these quantities, which enter the effective cross section (3.5), are needed to estimate the DVLF pair annihilations and, ultimately, their impact on the DM energy density.As we have done for the above-threshold states, we rely on the pNRQCD framework to obtain the relevant observables for the bound-state dynamics.A comment is in order.Boundstate decays are the counterpart of the local four-fermion operators projected onto bound states rather than scattering states (see detailed discussions in ref. [64]).Hence, the same arguments about the applicability of NRQCD and pNRQCD that involve the decay of heavy DVLF pairs into light SM states and/or scalar leptoquarks applies also here.
Bound-state formation:-For the model at hand, there exist two classes of boundstate formation processes, and corresponding decays, that depends on the DVLF pairs.One find color-singlet bound states, which originate from the combination of a DVLF particle and antiparticle.Moreover, by assembling two vector-like particles color-triplet bound states that appear together with a continuum spectrum of unbound pairs in a color sextet configuration.As for ψψ annihilations in a scattering state, see eq. (3.11), the counterpart for the negative-energy part of the spectrum is the decay of a spin-triplet color-triplet bound state (cfr.eq.(3.25)).As done in section, we do not distinguish explicitly between the color representation of ψ ψ pairs and their conjugate color pairs ψψ, and simply refer to them as color-triplet bound states and color-sextet scattering (or unbound) states.
The first ingredient for the estimation of bound-state effects on the DM energy density is the determination of the cross section for the two following processes where a color-singlet bound state B n [1] is formed from a color-octet scattering state (ψ ψ) p [8] via the emission of an ultrasoft gluon. 17The same hold for the second process that involves unbound pairs in a color-sextet (ψψ) p [6] and a bound state in a color triplet B n [3] .The subscripts indicate the color representation of the pairs, n stands for the collective discrete quantum numbers of a given bound state (n m), and the unbound scattering state is labeled with the momentum of the relative motion p = m ψ v rel /2.
The bound-state formation cross section can be computed at leading order from the imaginary part of the one-loop self energy in pNRQCD.This has been recently discussed 17 Bound states cannot form by the gluon emission from color-singlet or color-antitriplet scattering states because of SU(3) charge conservation, namely the two processes (ψ ψ) n and detailed in refs.[62,64] for abelian and non-abelian dark matter models (see ref. [67] for the case of soft scalar exchange and the corresponding pNREFT).We show the diagrams for the processes in eq.(3.13) in figure 5.In pNRQCD transitions among pairs are induced, at leading order, by chromoelectric-dipole vertices.The bound-state formation process (ψ ψ) p [8] → B n [1] + g has been computed formerly in the literature [29][30][31]64].We instead have to derive the corresponding process for transition from a color-sextet unbound state to a color-triplet bound state.We shall follow the procedure outlined in ref. [64].
Because we are interested in the bound-state formation process happening in the early universe, the gluon can be thermal, and, therefore, the computation needs to be performed in the thermal field theory version of pNRQCD.As long as the temperature scale is not larger than the inverse Bohr radius, one can rely on the in-vacuum derivation pNRQCD (see [58,68] for heavy-quarkonium and for dark matter [64]).In the so-obtained EFT, the dynamical scales are the ultra-soft scale, m ψ α 2 s , and the temperature scale.The main relevant aspect is that transitions among pairs are still described by the in-vacuum electric dipole transitions. 18 In order to compute the bound-state formation process (ψψ + g, we use pNRQCD for two heavy quarks [57], which well applies to the vector-like quarks of our model.The non-relativistic DVLF fields define a pair in the color space as follows where T (r, R, t) and Σ σ (r, R, t) are the bi-local fields of pNRQCD that accounts for the wave-function of the corresponding color configurations, r ≡ x 1 −x 2 is the distance between a fermion located at x 1 and an antifermion located at x 2 and R ≡ (x 1 + x 2 )/2 is the center of mass coordinate.Then = 1, 2, 3, σ = 1, ..., 6 and i, j = 1, 2, 3 and the tensors T ij and Σ σ ij can be found in ref. [57].The relevant electric-dipole interactions read [57] where T a are the SU(3) generators.Having clarified the relevant vertices, we move to the evaluation of the sextet self-energy in figure 5. We use the real-time Schwinger-Keldysh formalism [69,70].The real-time formalism leads to a doubling of the degrees of freedom called of type 1 and 2. The type 1 fields are the physical ones, namely those that appear in the initial and final states.Propagators are represented by 2 × 2 matrices, as they may 18 The multipole expansion holds for thermal gluons as long as the typical distance of the fermionantifermion or fermion-fermion pairs is smaller than 1/T .At large temperatures, T m ψ αs, the multipole expansion breaks down.In our numerical study, we solve the Boltzmann equation (3.1) starting from mχ/T = 10 down to smaller temperatures.For αs ≈ 0.1, the multiple expansion holds to a large extent for the whole temperature window, including the chemical freeze-out occurring for mχ/T ∼ 1/25.involve fields of both types.As for heavy non-relativistic particles at finite temperature, it has been shown in [58] that the 12 component of a heavy-field propagator vanishes in the heavy-mass limit, as a result, the physical heavy fields do not propagate into type 2 fields.Hence, the type 2 fermion-antifermion fields decouple and may be ignored in the heavymass limit, which reduces the relevant self-energies that we need to compute.In practice, it suffices to obtain the 11 component of the self-energy diagrams given in figure 5.As for the thermal gluon propagator, we shall adopt its form in Coulomb gauge [58,64].
The self-energy for the color-sextet field Σ σ (r, R, t) in the right panel of figure 5, reads19 where p 0 is the energy of the incoming pair and n B (x) = 1/(e x − 1) the Bose-Einstein distribution.In eq.(3.16) one can distinguish the in-vacuum and thermal contributions originating from the gluon propagator.The next steps are to extract the imaginary part of the self-energy in eq.(3.16), project the self-energy onto external scattering states and use the optical theorem in order to obtain the corresponding cross section Finally, we project onto intermediate color-antitriplet bound states, introduce a short-hand notation for the process (ψψ + g with simpler subscript and superscript for the cross section, and perform the color average of the cross section.Our result reads where we make explicit that one power of the strong coupling constant is evaluated at the ultrasoft scales µ us ≡ m ψ α 2 s , T , as dictated by the ultrasoft interaction in eq.(3.15) .Then, the energy difference between the incoming scattering state and outgoing bound state is at leading order where we used C a = (N + 1)/(2N ) in order to write the Coulombic energy levels in a compact way.As a relevant example, which we shall use in the numerical extraction of the DM energy density in section 3.3, we specify the general result in eq.(3.18) to the formation of the lowest-lying 1S bound state.In this case, only scattering states in the partial wave = 1 contribute, and the bound-state formation cross section reads (σ In the dipole matrix element, the natural renormalization scale of the coupling is the soft scale 20 .However, in order to avoid clutter, we dropped the corresponding scale dependence for the strong coupling in the bound-state formation expression (3.20).Details on the Coulombic wave functions and the general expression for the electric-dipole matrix elements for sextet-to-antitriplet transitions are given in appendix B.
In order to compare with the bound-state formation process (ψ ψ) p + g, that is also needed for the determination of the DM energy density in section 3.3, we provide the cross section for the formation of the color-singlet ground state, which reads [29,31] (see refs.[30,64] for a derivation in pNRQCD) (3.21)The corresponding energy difference between the incoming color-octet scattering state and the color-singlet bound state is In figure 6 we show the thermally averaged bound-state formation cross section for the 20 We do not distinguish the soft scale between antitriplet bound states and unbound sextets.color-singlet and color-triplet 1S bound states divided by πα s (µ h ) 2 /m 2 ψ . 21The thermal average is performed in the standard way, see e.g.[19,64], with Maxwell-Boltzmann distribution of the incoming DVLF pair.On the one hand, one may see how the two different bound-state formation processes are comparable at typical freeze-out temperatures, with the cross section for (ψ ψ) p [8] → B n [1] + g being marginally larger than (ψψ) n [6] → B n [3] + g.On the other hand, at smaller temperatures, the bound-state formation for the color-triplet 1S state is larger and peaks at later time with respect to the color-singlet ground state.This latter aspect is due to a smaller absolute value of the binding energy for the triplet, as for the latter aspect, a qualitative similar behaviour is found when comparing the ground state with excited states [30,51,64]).
Dissociation and decays:-Once bound state form, it can either decay or get dissociated by thermal particles of the early universe plasma.In this work, we consider the gluodissociation process [68,71,72], namely . Whenever thermal gluons in the early universe plasma have sufficient energy, they can break the bound state and turn it into an above-threshold scattering state.The corresponding rate is a thermal 21 The normalization factor is just needed to plot a dimensionless quantity, Sbsf ≡ (σ bsf v rel )/(παs(µ h ) 2 /m 2 ψ ).One may consider different renormalization scale for the strong coupling in the normalization factor, e.g.παs(µs) 2 /m 2 ψ , which would make the curves lower.
width, or dissociation width, of a bound state.The efficiency of the conversion of bound states into its decay product depend on the interplay of the dissociation and decay width, as displayed in the effective cross section in eq.(3.5).The bound-state dissociation can be obtained in two ways.One powerful argument is that, whenever the ionization equilibrium is maintained, the bound-state dissociation and bound-state formation cross section are related via the Milne relation [19,29] (one can find a recent derivation for it in ref. [31]).More specifically, the Milne relation links the bound-state formation with the ionization cross section.The latter is used to obtain the bound-state dissociation width through a convolution integral with the incoming thermal gauge boson momentum, that may break the bound state if sufficiently energetic.We write the dissociation width for a generic color configuration [R] = [1], [3], as follows where |k| is the energy of the gluon, [R] stands for the representation of the bound state either in a color singlet or a color antitriplet, g g are the gluon degrees of freedom and σ n ion,[R] (|k|) is the ionization cross section, which is related to the bound state formation cross section via the Milne relation Here g B [R] stands for the degrees of a bound state, whereas [R ] for the unbound scattering states in a color-octet or color-sextet represenation.Alternatively, and as a non-trivial check, one can derive the bounds-state dissociation from the imaginary part of the bound-state in pNRQCD.This has been shown in refs.[64,67] for vector as well as for scalar force mediators in the context of dark matter, and earlier in refs.[58,68] for heavy quarkonium phenomenology.
We do not include the additional dissociation mechanism as induced by 2 → 2 scatterings with the in-medium constituents.This is known as inelastic parton scattering [73][74][75][76] in heavy quarkonium literature.The counterpart for cosmological applications to dark matter freeze-out was considered in refs.[50,63,77] in the screening regime, namely for temperature larger than the typical Bohr radius T m ψ α s .Here, a non-trivial interplay with another thermal scale, a Debye mass for the gluons, is established and the extraction of the relevant cross sections and widths of the pairs becomes rather challenging for a broad temperature range.The bound-state formation process has been computed at fixed order, without a resummation of collective plasma effects that generates a Debye mass for the gluons, in ref. [78] (the corresponding cross section for an abelian dark matter model is given in ref. [79]).A careful investigation of the Debye mass scale within the framework of non-relativistic effective field theories, in particular its role in the bound-state formation via gauge boson emission, is still ongoing and we do not account for it in our work (however see [30] for an exemplary implementation of these effects in a dark matter model with colored coannihilators).Hence, the DM energy density as derived in section 3.3 has to be understood as upper bound.
The last ingredient is the decay width of the bound states.The color and spin-averaged bound state decay widths for DVLFs pairs, at leading order in the coupling and in the nonrelativistic expansion, read as follows ,ann stands for decay width of nS color-singlet bound states, which receive contributions from decays into gluon and leptoquark pairs, whereas Γ n [3],ann is the color-triplet nS decay width, which encompasses only decay into leptoquark pairs.We have explicitly indicated the scale for α s at the hard scale, µ h = 2m ψ , and the soft scale µ s , which originates from the wave function.
Finally, in the right panel of figure 6, we show the effective bound state formation cross section, namely the second term in the right-hand side of eq.(3.5), once again normalized by πα s (µ h ) 2 /m 2 ψ in order to display a dimensionless quantity.Here, the dissociation widths as well as the decay widths enter.Further parameters are specified in the plot label.One can see how the bound-state formation of a color-triplet 1S state (orange-dashed curve) gives a sizeable contribution to the total bound-state formation (red dash-dotted line) for y = 1.0.We have checked that the (ψ ψ) p [8] → B n [1] + g is largely dominant for y 0.5.

Numerical results for the DM energy density
In this section we solve the effective Boltzmann equation (3.1) with the relevant cross sections and decay widths that have been discussed in sections 3.2.1 and 3.2.2.The main scope is to address the impact of near-threshold effects on the DM energy density in the coannihilation regime, namely for small mass splittings between the dark fermion χ and the DVLF.We recall at this stage that the Sommerfeld and bound-state effects play a role for the DVLF pair annihilations, namely when the incoming states are ψ ψ, ψψ and ψ ψ.In figure 7 the dark matter energy density is given as a function of the dark matter mass m χ for the two options, Majorana and Dirac dark matter, respectively in the left and right panel.We take the leptoquark mass to be m φ = 1.5 TeV, the relative mass splitting is ∆m/m χ = 10 −3 and the portal coupling is y = 1.5.In this section, we fix the leptoquark-to-SM couplings y L = y R = 1 (we elaborate more on varying these couplings in section 4.2).The orange-dotted, purple-dashed and brown lines correspond to the energy density as extracted with free annihilation cross section, the Sommerfeld-only corrected cross section and with the inclusion of both Sommerfeld effects and bound-state formation respectively.More specifically, the solid-brown curve accounts for the lowest lying 1S bound state, whereas the brown-dashed line comprises the effect of the 2S state as well, in the no-transition limit.In the Majorana case, both bound states in a color singlet and color antitriplet contribute.
In addition to the dip at m χ ≈ m φ /2, which captures the resonant enhancement of the coannihilation process χψ → q c ¯ , one finds a further dark matter mass range where the energy density is locally decreased.This is due to the opening of additional annihilation channels for the dark fermion χ and DVLF into leptoquark pairs.In the Majorana dark matter scenario, the Sommerfeld effect as well as formation and decays of bound states make this feature more prominent.As a general trend, below the leptoquark mass threshold, the overall Sommerfeld corrections have practically no impact (a small effect can be seen at the resonant window m χ ≈ m φ /2).One can trace this back to a competing enhancement and suppression of the color-singlet and color-octet contributions to the annihilations into light Standard Model QCD states, cfr.eq.(3.7), that makes the Sommerfeld-corrected cross section slightly smaller than the free cross section at the freeze-out.However, the situation is different above the leptoquark mass, where DVLF annihilations experience an overall enhancing effect (accordingly the purple-dashed line is below the orange-dotted line because of a larger cross section that results in a smaller DM energy density).For the Majorana option, the enhancement of the cross section is more important because of ψψ → φφ and its conjugate process, whose leading contribution to the annihilation cross section originates from an attractive color-triplet channel, see eq. (3.11).As for the Dirac case, only ψ ψ → φφ † can occur, for which competing color-singlet and color-octet effects make the cross section only slightly larger than the free cross section, see eq. (3.9).The bound-state effects have a rather different behaviour with respect to the Sommerfeldonly scenario, see solid-brown lines in figure 7. Bound-state formation is effectively active below the leptoquark mass threshold due to formation of color-singlet bound states and their decays into light SM quarks and gluons.Above the leptoquark mass, the formation and decays of color-antitriplet bound states also contribute in depleting the dark matter, since the corresponding decays becomes kinematically allowed. 22For m ψ larger than the leptoquark mass, color singlet bound-state can also decay into four-body SM states through unstable φφ † pairs.Since bound-state effects are efficiently annihilating the DVLF pairs, and hence the DM fermion in the coannihilation regime, the DM energy density is systematically below the free-annihilation scenario for the entire dark-matter mass range.As one may see from the comparison of the left and right plots in figure 7, bound-state effects are more important for the Majorana case, because of the additional bound-state formation processes for (ψψ) and ( ψ ψ) pairs and corresponding bound-state decays.
For the choice of the parameters as given in figure 7, and accounting for 2S excited state in the non-transition limit, the dark matter mass that is consistent with the observed energy density shift from m χ = 3.4 TeV (m χ = 3.1 TeV) to m χ = 5.6 TeV (m χ = 3.9 TeV) for the Majorana (Dirac) case.Having clarified the impact of the Sommerfeld-only corrected cross section, in the following we present the numerical results by accounting for both manifestation of non-perturbative effects, i.e.Sommerfeld and bound states.Changing the leptoquark mass does not affect the qualitative behavior of the DM energy density curves displayed in figure 7.
In the next example, we consider the DM energy density as a function of the leptoquark mass.In figure 8, left panel, we notice how the non-perturbative effects provide a much wider range for m φ that is compatible with Ω DM h 2 ≤ 0.1200 ± 0.001, so that we do not overclose the universe up to m φ = 5.6 TeV.On the contrary, if one estimates the DM energy density without Sommerfeld and bound-state effects, leptoquark masses 2.4 TeV ≤ m φ ≤ 4.5 TeV are excluded by the Planck collaboration measurement, and one has to rely on the resonant dip that allows for a viable mass range 4.6 TeV ≤ m φ ≤ 5.6 TeV.For the same choice of parameters, the Dirac case displays differences with the Majorana dark matter scenario.The overall annihilation cross section is smaller and the curves shift at a higher DM energy density.The less prominent non-perturbative effects make the orange-dotted line and brown curves closer for m χ and m ψ larger than m φ .Moreover, the leptoquark mass window compatible with Ω DM h 2 ≤ 0.1200 ± 0.001 is m φ ≤ 1.0 TeV for the free case, whereas it is extended to m φ ≤ 1.5 when Sommerfeld and bound-state effects are included.Then, a second mass region is available at around m φ ≈ 2m χ , as the resonant enhancement is sufficient to reduce the DM energy density below the observed value.
Finally we aim to explore the parameter space of the model which is compatible with the observed dark matter energy density.The model contains three mass parameters (m χ , m ψ , m φ ) and three new couplings (y, y L , y R ).We focus on (m ψ , m φ ) contours for different values of the relative mass splitting ∆m/m χ and the portal coupling y (we remind the reader we fix the leptoquark-to-SM couplings y L = y R = 1).Our choice to consider the (m ψ , m φ ) mass plane is motivated by the present collider limits that are largely applicable to the colored states of the model, see discussion in section 4.3.
In figure 9 we provide the curves that reproduce the observed energy density for Majorana (upper row) and Dirac (lower row) dark matter.Let us discuss the first scenario.We fix the relative mass splitting to ∆m/m χ = 0.05 in the left panel and to ∆m/m χ = 0.01 in the right plot.We select three benchmark values for the portal coupling to be y = 0.1, 0.5, 1.0.The qualitative difference between the left and right panels of figure 9 can be explained as follow.The mass splitting ∆m/m χ = 0.05 makes the role of ψψ, ψ ψ and ψ ψ annihilations less important than the case with ∆m/m χ = 0.01.Accordingly, in order to reproduce the same energy density, higher dark matter mass, and hence m ψ , are required.As a result the curves are shifted towards larger DVLF masses, irrespective of the portal coupling y.Another aspect worth explaining is the behaviour with varying y.The effect is well visible in both panels of figure 9 (upper row).By increasing the value of y, it corresponds to a larger cross section of the resonantly-enhanced annihilation ψχ → q c ¯ /Q c L, with the leptoquark in the s-channel.For m χ ≈ m ψ ≈ m φ /2, the enhancement demands large dark matter and DVLF masses in order to reproduce the observed energy density.The effects gradually fades away for decreasing y.For large m φ masses, the contours merge into straight vertical lines independent of y, which signals that the energy density is determined by DVLF annihilations into SM quarks and gluons.The gray shaded areas implements the relevant collider exclusion limits, m φ , m ψ ≤ 1.5 TeV.For ∆m/m χ = 0.05, a good portion of the cosmologically favoured parameter is probed and ruled out for small y's.The surviving regions are those along the resonant condition, m ψ > m φ , and for large Yukawa-portal coupling y = 1.For the smaller relative mass splitting ∆m/m χ 0.01, the main constraint comes from the leptoquark exclusion limit, since the required DVLF masses are m ψ 2.2 TeV and are out of the reach of current collider limits.
The corresponding parameter space for the Dirac dark matter option is displayed in figure 9 (lower row).One find the main same qualitative features as for the Majorana case in figure 9.However, since the total annihailtion cross is smaller in the Dirac case, the parameter space compatible with the observed DM energy density shifts to smaller m ψ ≈ m χ masses.As a result, for the larger relative mass splitting ∆m/m χ = 0.05, only the case y = 1.0 remains still viable, and only along the resonant region.For the smaller splitting, the stronger coannihilations also make the DVLF masses out of the present collider limit, though in a less severe way with respect to the Majorana option.

Correlated observables, dark matter and collider phenomenologies
In this section, we discuss the interplay among different observables as well as summarize the collider constraints on the masses and the relevant Yukawa couplings.

Dark matter direct and indirect detection prospects
The DM fermion χ can interact with the nucleon constituents only via loop processes, where the DVLQ and the leptoquark run in the loops.Following ref. [80], we calculate the spinindependent cross-section σ SI for the Majorana (Dirac) DM option at one-loop, and we find that the typical cross section are beyond the current, and most likely future, sensitivities.For example, for the benchmark point satisfying the correct relic abundance as given in figure 7, namely m χ = 5 (3.5)TeV, ∆m/m χ = 0.001, m φ = 1.5 TeV, the corresponding cross-section is σ SI = 2.6 (4.7) × 10 −51 cm 2 , which is not only several orders of magnitude below the current experimental limit, 1.45 × 10 −45 (9 × 10 −46 ) cm 2 [81], but also below the neutrino coherent scattering floor, σ νN ∼ 10 −48 cm 2 [82].Therefore, both Majorana and Dirac dark matter candidates are beyond the reach of direct detection experiments in the foreseeable future.
An additional potential signal is given by the present-day annihilations of DM particles from astrophysical sources, that can be searched with indirect detection strategies.In our case, Majorana or Dirac DM can leave an imprint via annihilation into leptoquark pairs and their subsequent decays into the SM quarks and leptons at the tree-level, namely χχ(χ χ) → φφ * → q q , as well as via loop-induced processes, i.e. χχ(χ χ) → γγ, gg, Z Z, γ Z. Here, we simply give an estimation of the cross section for dark matter being heavier than the leptoquark.In this case, the tree-level DM annihilation into four-body final states via the decay of unstable leptoquarks is not phase-space suppressed, at variance with the offshell region m φ > m χ , whereas the loop-processes remain suppressed.For Majorana DM, the corresponding pair annihilation into leptoquarks is p-wave suppressed (cfr.eq.(A.6)) and, therefore, the present-day annihilation rate of Majorana DM will be several orders of magnitude smaller than the case of Dirac DM.We then focus on the latter option in the following.In figure 10 (left), we present the present-day Dirac DM annihilation cross-section into on-shell leptoquarks for a fixed value m φ = 1.5 TeV and for (m χ , y) pairs that give the correct relic abundance.Two relative mass splittings are considered.Given that DM fermion masses lie in the TeV range, very high energy (VHE) gamma-rays can be expected from the energetic SM final state particles as produced from DM annihilation in typical DM-rich environments like the Galactic Center (GC) and Dwarf Spheroidal (dSph) galaxies.For DM in the mass range of 1-10 TeV, the current combined limits from VHE gamma-rays searches at 20 dSph galaxies by the Fermi-LAT, HAWC, H.E.S.S., MAGIC, and VERITAS experiments [83] on the DM annihilation are in the range 1.2×10 −25 −3.3×10 −24 cm 3 /s for the bb final state, and 4.4 × 10 −25 − 4.3 × 10 −24 cm 3 /s for the τ + τ − final state, respectively.Moreover, the future sensitivities of the upcoming Cherenkov Telescope Array (CTA) on the bb and τ + τ − final states for the Galactic Center with Einasto DM profile and an observation time of 525 hr [85], are (1.4 − 2) × 10 −26 cm 3 /s and (3 − 9) × 10 −26 cm 3 /s, respectively, for 1-10 TeV DM mass.
Despite the annihilation cross-section curves cover the same order of magnitude of the current combined limits or future sensitivities, we stress that the experimental limits are derived for two-body final states.In our case, DM annihilating to leptoquark pairs produces in turn four-body final states consisting of two quarks and two leptons (accordingly we refrain from superimposing the experimental limits with the cross section curves of the model).As we can see from figure 10 (right), for a banchmark value m χ = 3.5 TeV, the end-point energy and photon spectra differ for the four-body final states, t τ t τ , bνbν, and t τ b ν + h.c compared to the two-body final states, bb and τ + τ − .Therefore, a dedicated statistical analysis is required to quantitatively derive the limit from the current experiments or the sensitivity of the future experiments on the Dirac DM annihilation of the model, and assess the indirect detection prospects.

Muon anomalous magnetic moment
Lepton flavor universality is not a fundamental property of nature.New physics can, in principle, couple more strongly to a specific fermion generation.In fact, there is a longstanding tension in the muon anomalous magnetic moment, a µ = (g − 2) µ /2.This discrepancy was measured at the E821 experiment [86] in 2006, which was recently confirmed by the E989 experiment [87].The combined result yields a 4.2σ discrepancy with the SM prediction, hinting towards physics beyond the Standard Model that violates lepton flavor universality.
Interestingly, the example model that we have discussed in the previous sections contains a scalar LQ, φ ∼ (3, 1, 1/3), and it can address the tension in the muon anomalous magnetic moment (for an incomplete reference list, see, for example, refs.[88][89][90][91]).Since φ couples to both the left-and the right-handed up-type quarks, it is possible to have a chiral enhancement in the loop to provide adequate new physics contributions to (g − 2) µ .
To compute the (g −2) µ , we work in the up-type quark mass diagonal basis (for details, see ref. [90] and references therein), where the CKM matrix is associated with the downtype quarks.In this basis, the Yukawa couplings of the LQ, cfr.eq.(2.5), take the following form: where V represents the CKM matrix; for its entries, we use the PDG values [92].With the couplings as given above, additional contributions to the (g − 2) µ are generated at the one-loop level, which can be expressed as follows [90]: here we have defined, x q = m 2 q /m 2 φ , and the index q runs over u, c, t-quarks.A sufficient new physic contribution to the (g − 2) µ can only be provided if a top-quark or charm-quark mass flip occurs inside the loop.Due to the very small mass, the contribution from the upquark can be fully neglected.This is why we examine two separate scenarios, (i) one with y L,R tµ = 0 (top-quark mass flip) and (ii) another with y L,R cµ = 0 (charm-quark mass flip).Moreover, in eq. ( 4.3), the first term corresponds to the chirality-flipping contribution, hence it dominates over the second term, which can be safely neglected.Then, the (g − 2) µ becomes approximately proportional to the ratio y L qµ y R qµ m q /m 2 φ .Since m t /m c ∼ 136, with ) For simplicity, here we quoted the most stringent collider bound (gray-shaded area), which can be relaxed depending on the model/analysis details, see discussion in the text.The purple band gives the observed DM energy density for m χ = 2.5 TeV and ∆m/m χ = 10 −2 .
, the experimentally measured deviation in the muon anomalous magnetic moment can be incorporated for m φ 5 TeV and m φ 60 TeV for the charm-and top-quark scenarios, respectively.Note that due to the appearance of the CKM matrix in eq.(4.2), all three generations of down-type quark couple to the (muon) neutrino.Depending on the scenario, additional interactions of these types may lead to uncontrollable flavor-violating processes (for example, an additional interaction with the electron may mediate dangerous µ → eγ processes [90]).Owing to the requirement of small couplings for TeV scale LQ in the topquark mediated scenario, flavor-violating processes (for example, B → K ( * ) νν) are well under control.This scenario with top-quark mass flip is illustrated in figure 11.
The charm-quark mediated case, however, suffers from large flavor-violating processes.This is due to the requirement of somewhat larger couplings.The scalar LQ couples to strange-and down-quark with almost the same strength, where the latter coupling is Cabbibo suppressed.As a result, the leptoquark mediates kaon decays of the type K + → π + νν, which rules out a large part of the parameter space as shown in figure 12 (green-shaded area).The experimental result from NA62 [93] that corresponds to BR K + → π + νν < 1.85 × 10 −10 puts strong constraint only on y L cµ for our scenario, which we compute following Ref.[94].Therefore, in order to get the correct (g − 2) µ value as observed in the experiments, a large y R cµ is typically required.Interestingly, such a large value of y R cµ is also highly constrained [95,96] from non-resonant dilepton searches at the LHC [97,98], to be discussed below, see section 4.3.Once these two constraints are imposed, a tiny portion of the parameter space remains consistent with the experimentally observed (g − 2) µ only at the 2σ C.L. as depicted in figure 12. ).The green-shaded area corresponds to the experimental result from NA62 [93], whereas the gray-shaded area to non-resonant dilepton searches at the LHC [97,98].The purple band gives the observed DM energy density.
In order to highlight the interplay between the explanation of the observed values (g − 2) µ and dark matter relic density, we superimpose the curves that reproduce Ω DM h 2 = 0.1200 ± 0.001 in figures 11 and 12.To this aim, we consider the Majorana dark matter option.As regards the charm case, the quite small available parameter space demands a careful choice of the DM mass once the portal coupling y and the mass splitting are fixed; we choose y = 1 and ∆m/m χ = 10 −2 .For the right-handed Yukawa value y R cµ = 1.40 (y R cµ = 1.96), we find the corresponding DM mass to be m χ = 2.84 (2.64) TeV in order to lie in the still viable window.Changing the DM mass more than about 5% is sufficient for loosing the interplay (one could however tune again the mass splitting and y).The top quark scenario is much less restrictive in this respect and we show an exemplary case, which lies beyond the most stringent LHC limit but it is still in the TeV range of the leptoquark mass.Here, the DM is m χ = 2.5 TeV, whereas we keep the same values for y and ∆m/m χ as in the charm case.
Before concluding the discussion about the muon anomalous magnetic moment, we point out that recent lattice determinations [99][100][101] of the hadronic vacuum polarization give a SM prediction that agrees with the experimental result, however, it is in tension with the previous calculations based on dispersive methods [99].Forthcoming experiments will be able to shed light on this unresolved issue.

Collider implications
LHC constraints:-Since the DM in our framework is a SM singlet, it cannot be directly produced at the LHC.However, the corresponding productions of the LQ and the DVLF are unsuppressed since they carry color charge.Relevant bounds on these masses from LHC searches are discussed in the following.
Leptoquarks can be pair produced at LHC via gluon-fusion pp → φφ † [102,103].After production, each of these LQs would decay into a quark and a lepton.Several dedicated searches for LQ pairs have been carried out by ATLAS and CMS Collaborations for different final states with (pp → qqνν) or without (pp → qq ) neutrinos.The LHC limits on LQ mass depend on the branching ratios of various modes.For the top-quark mass flip solution presented above, the strongest (weakest) constraints arise if the value of y R tµ (y L tµ ) is somewhat larger than y L tµ (y R tµ ).In this case, branching ratio to pp → ttµµ is about unity (a half) and LHC excludes LQ masses below 1.5 (1.3) TeV [96,104].
Note that when kinematically allowed, in addition to the φ → q , there is an additional decay mode, namely, φ → ψχ (see eq. (3.4)).As a result, the branching ratio to q gets modified and the bounds quoted above will be relaxed (for details, see Refs.[35,36]).For pair produced LQs, a dedicated search has been performed at the LHC, where one of the LQs directly decays to jet and a muon (φ → jµ) and the other LQ decays to φ → jµ + / E T with low-p T SM fermions via cascade decays (φ → ψχ → φ * χχ → jµχχ).The analysis strategy is based on the search of a peak in the LQ invariant mass m jµ distribution from the highest p T muon and jet in an event, with the requirement of significant missing transverse momentum due to the DM particles in the final state.Since no signals above the SM background is observed, from this search, LHC rules out dark matter masses up to 600 GeV for LQ masses of order O(1.5)TeV [105].
On the other hand, for the charm-quark mass flip solution, the most relevant LHC bound comes from the indirect high-p T searches [106].We are interested in the direct constraint on the coupling versus mass plane arising from the non-resonant dimuon searches at the LHC (pp → µµ).As discussed above, we require somewhat large values of y R cµ to address the tension in the muon magnetic moment.For a LQ of mass 1 TeV, non-resonant dimuon search rules out couplings of order unity, i.e., y R cµ ≤ 1.1 [95] must be satisfied.Since the bound on the coupling depends on the mass of the LQ, this functional dependence is presented in figure 12 with varying m φ for the two different coupling choices (gray shaded area).
DVLFs are also pair-produced at the LHC via gluon-fusion pp → ψψ.Subsequently, each DVLF decays to ψ → χ φ → χ q leading to large MET.Processes of these types have been searched for at the LHC that put strong bounds on the lower limit of squark masses, especially for stop and sbottom.Depending on the exact LQ coupling and mass, as well as the value of ∆m = m ψ − m χ , the bound is in between m ψ 700 GeV and m ψ 1300 GeV (for details, see Ref. [36]).This analysis, however, is not applicable for very small mass splitting.In such a compressed scenario, missing energy searches lose sensitivity, and conventional searches are no longer applicable.This happens typically for mass splittings ∆m < 5 GeV [107].In fact, if the mass splitting is very small, the DVLF becomes effectively long-lived as a result of a highly off-shell LQ.The phenomenology of DVLFs is entirely different from the one in the standard searches of these particles.These quasi-stable heavy quarks, namely, the R-hadrons [108] interact hadronically as they move through the detector after being produced at the LHC.The recent analysis of the ATLAS collaboration puts bounds on the mass of long-lived supersymmetric R-hadrons (squarks and gluinos), which for a state with electromagnetic charge ±1/3 (sbottom) is 1250 GeV [109].This search is quite model-independent and has been adapted to the case of vectorlike fermion, see Ref. [110], which finds a lower bound of m ψ 1500 GeV.
It is important to point out that the exact LHC limits depend on the details of the multidimensional parameter space, which is beyond the scope of this work.
Muon collider probes:-As discussed above, the observed large tension in the muon anomalous magnetic moment is an indication that the new physics couples strongly with the muon and not to the other lepton generations.Consequently, muon colliders are the perfect machines to test such muon-philic new physics scenarios [111].As outlined above, the LQ can reside in the multi-TeV range, which is beyond the reach of LHC, and yet provide the required new physics contribution to reproduce muon g − 2 and play an important role in determining the dark matter relic abundance.In such a scenario, by integrating out the heavy field in obtaining the effective field theory, one gets the scattering process µµ → cc/tt depending on the charm-philic/top-philic nature of the LQ.Then, a probe of ∆a µ is obtained via computing the number of events and requiring a statistically significant deviation from the SM µµ → cc/tt background.By performing extensive analysis for the relevant semi-leptonic operator involving charm-quark (top-quark), ref. [111] showed that muon-philic scenario can be probed already at √ s = 1 TeV, while the top-philic case can be probed at √ s = 10 TeV.Such a high-energy determination of ∆a µ is a unique feature of muon colliders.

Conclusions
In this work, we considered a class of dark matter models where the DM candidate does not interact directly with Standard Model particles.The sole interaction of the dark matter, a SM singlet fermion, is through a Yukawa coupling with a scalar mediator and a Dirac fermion.The latter is assumed to carry a dark charge and to be heavier than the DM particle to ensure the DM candidate's stability.In order to trigger DM annihilations in the early universe via thermal freeze-out, both the mediator and the fermionic partner carry some charges under the SM gauge group.We focused on the case of QCD colored states and the scalar mediator being a leptoquark.This setup provides interesting phenomenological consequences despite the absence of direct interaction between the dark matter and the visible sector-most notably, collider implications and a connection with the anomalous magnetic moment of the muon.
The colored DVLFs play a crucial role for the annihilations of particles of the dark sector, that determine the relic energy density.Indeed, in order to avoid overclosing the universe, coannihilations of nearly degenerated colored partners have been shown to be a necessary ingredient for the model at hand.In this regime, it is important to scrutinize relevant non-perturbative effects.One of the main objectives of this paper is to assess such effects for a more reliable estimation of the dark matter energy density.
Non-relativistic heavy DVLF pairs are affected by repeated soft-gluon exchange in two ways.First, above-threshold scattering states experience Sommerfeld effects.We show that they play a rather marginal role in the case of ψ ψ annihilations because of a competing enhancement and suppression in the attractive color-singlet and repulsive coloroctet channels, that contribute at the same order in the velocity expansion (see eqs. (3.7) and (3.9)).On the contrary, a more prominent effect is found for ψψ and ψ ψ annihilations because of an enhancing color-triplet Sommerfeld factor.The corresponding cross section is not diminished by a suppression factor from the color-sextet repulsive channel, which only appears at higher order in the velocity expansion.This latter situation only applies to the Majorana dark matter option.Second, repeated gluon-exchange in the attractive color-singlet and triplet channels, may sustain bound states.During the freeze-out in the early universe, bound-state formation for DVLF pairs and their subsequent decays into SM particles, works as an additional channel to effectively deplete dark sector particles.We take into account the bound-state formation process via gluon radiation.For the Majorana fermion scenario, we have obtained the bound-state formation cross section for the process (ψψ) n [6] → B n [3] +g in the framework of pNRQCD, and computed the corresponding electricdipole matrix elements in full generality (see appendix B).Our result can be also be useful for other simplified models that feature real scalar dark matter coannihilating with vectorlike colored fermions, e.g.[112,113].
We have assessed the impact of non-perturbative effects depending on the nature of the dark matter fermion, which is rarely pursued in the literature.As a general observation, Sommerfeld and bound-state effects are more relevant for the Majorana option.This is due to additional DVLF pair annihilation channels, namely ψψ → φφ and the complex conjugate process.For this scenario, bound-state formation and decays from both color-singlet and color-triplet pairs boost the annihilations of the dark sector particles.The inclusion of non-perturbative effects, especially bound-state formation, has a sizeable impact on the dark matter mass that is compatible with the observed energy density.For the smallest mass splitting considered in this work, we find that m χ is shifted from 3.4 TeV to 5.6 TeV (see figure 7).Moreover, as shown in figure 8, non-perturbative effects open new mass regions for the leptoquark, which would be deemed excluded otherwise.Our findings motivate further investigations and a more comprehensive inclusion of bound-state effects for the models' class of this work (such as the complementary bound-state formation process via 2 → 2 scatterings with light plasma constituents and more excited states).
Despite we have considered the freeze-out option in this work, some comments can be made on the conversion-driven freeze-out.Here, much smaller portal coupling y are needed, which makes the Sommerfeld enhancement due to the attractive triplet channel, as well as the corresponding bound-state effects practically irrelevant.However, bound-state effects from (ψψ) n [8] → B n [1] + g, which are independent of the Yukawa coupling, can be relevant for estimating the thermal abundance of the DVLF partner both in the Dirac and Majorana option.We leave their inclusion in the conversion-driven freeze-out for future investigation on the subject.
Although dark matter direct and indirect detection is challenging in this setup, because of no coupling between the DM with the SM particles, the example model that we studied still has important phenomenological consequences.The mediator chosen is a scalar leptoquark, which can be directly searched for at colliders.For masses of order TeV, LHC already puts strong constraints on the LQ couplings.Interestingly, in addition to acting as a mediator between the visible and dark sectors, the leptoquark can also address the longstanding tension in the muon's anomalous magnetic moment.In such a scenario, leptoquark masses, even up to 100 TeV, can be probed in future muon colliders.Since the coannihilation partner must have identical quantum numbers as the mediator within this framework, it can also be efficiently produced in LHC and may leave exciting signatures.In view of the fact that coannihilations play a crucial role in achieving a dark matter abundance compatible with observations, which requires nearly degenerate states, detecting the coannihilation partner at LHC will already provide information about the mass of the dark matter.
where T is the temperature, K 1 (x) and K 2 (x) are the modified Bessel function of the first and second kind, λ(x, y, z) = (x − y − z) 2 − 4xyz and We remark that the analytical 2 → 2 cross sections have been checked against the model implementation in MadGraph [48].In the main body of the paper, in eq.(3.3), we have abbreviated the cross sections by indicating only the incoming states σ ab v rel .

χχ annihilations
The pair annihilation cross section of Majorana dark matter fermions into leptoquark φφ † pairs is where N is the number of colors, and we have defined the following auxiliary quantities in order to write the cross section more compactly The pair annihilation cross section of Dirac dark matter fermions into φφ † pairs reads instead Upon expanding in the non-relativistic velocity, which is a good approximation sufficiently away from the opening of the threshold, the cross section times the relative velocity reads for the Majorana case, whereas for the Dirac case we find We have kept here the corresponding leading terms in the velocity expansion.

χψ annihilations
In this case there are two class of processes: (i) annihilation processes into lepton and quark pairs via a s-channel leptoquark exchange; (ii) annihilations into a gluon and a leptoquark.In the latter case, the unstable leptoquark decays in turn into a lepton and a quark.The coannihilations into a right-handed lepton and quark reads where the factor of 2 simply originates from the SU(2) multiplicity.
For process χψ → φg the analytical expression of the cross section turns out to be quite lengthy.We list the squared amplitude of the s and t-channels, as well as the interference term.One can easily obtain the cross section by incorporating the flux factor, namely 2s λ(1, m 2 χ /s, m 2 ψ /s), the two-body final state (for a massive letpoquark and a massless gluon) and performing the spin and color averages.The s-channel squared amplitude is Particle-antiparticle annihilation of the DVLF are divided in two classes.On the one hand, there are 2 → 2 annihilation processes directly into light SM states, namely gluons and quarks.On the other hand, annihilation into leptoquark pairs are also viable.The latter induce a four-body final state which is relevant above the leptoquark threshold.For the first class, the velocity expansion works fine and one can readily extract the cross section from the matching coefficients of NRQCD [54] (see results in the body of the paper, cfr.eq.(3.7)).
As for the annihilation into leptoquark pairs, without performing the velocity expansion, we obtain (A.17) The result in eq.(A.15) enters the total cross section of both scenarios, namely a dark matter Dirac or Majorana fermion.Away from the leptoquark mass threshold, the expansion of the cross section in eq.(A.15) gives the result in eq.(3.9), when decomposed in the corresponding color singlet and color octet contributions.

ψψ and ψ ψ annihilations
When the dark matter fermion is assumed to be Majorana, there are additional annihilation channels for particle-particle (ψψ) and antiparticle-antiparticle ( ψ ψ) DVLF pairs.In this case, a t and u-channel mediated by the exchange of a Majorana fermion are possible.The cross section for the process reads where the auxiliary coefficeints E and F can be read off eqs.(A.16) and (A.17), and σ ψ ψ→φ † φ † (s) = σ ψψ→φφ (s).

B Electric-dipole matrix element
In this section we provide the derivation of the electric dipole matrix element for the transition between color-sextet scattering states and color-antitriplet bound states.We give the result for a generic bound state and by choosing p along the z-direction.We shall derive a general expression using the notation and decomposition of the scattering and bound-state wave functions following the derivation for the octet-singlet electric dipole (see ref. [64]), which was in turn based on refs.[115,116].The necessary ingredients are the wavefunctions of the Coulombic scattering and bound states.The Coulomb wavefunction for a DVLF scattering state of positive energy p 2 /M reads, when expanded into partial waves Ψ p (r) = r|p [6] as Ψ p [6] (r) =  (B.9) In the dipole matrix element the natural renormalization scale of the coupling α s is µ s , which is of the order of the soft scale.We do not distinguish the soft scale between singlet and triplet bound states.

Figure 1 .
Figure 1.Annihilation and coannihilation channels of the dark sector particles (from left to right).Dark matter annihilation χχ → SM SM (χχ → SM SM SM SM) takes place only (already) at the loop-level (tree-level).For the coannihilation partner, both the ψψ → SM SM and ψψ → SM SM SM SM occur at the tree-level.Moreover, coannihilation channels χψ into visible sector particles lead to SM SM and SM SM SM final states at the tree-level.Processes relevant for dark matter direct detection (for example, from bottom to top in the leftmost diagram for χ) only happen at loop-level.

Figure 2 .
Figure 2. Left: Representative diagram for DVLFs pair annihilation into light SM states; gluons are depicted with curly lines.Middle and right: DVLFs pair annihilations into leptoquarks; leptoquarks are depicted with dashed lines, whereas the solid double arrow stand for the DM fermion χ (forwardbackward arrow for χχ Majorana-type field contractions).Repeated gluon exchange are shown for the incoming DVLFs, both for ψ ψ and ψψ annihilations.

Figure 3 .
Figure 3. a) A representative diagram for the matching procedure that leads to four-fermion local operators in NRQCD.The imaginary part of the loop diagrams encode the annihiatlions of DVLFs into gluons with typical momenta |k f | ≈ m ψ .b) The box diagram captures the annihilations of DVLF into massive leptoquarks.For small mass splittings, the typical momenta of the final-state leptoquarks is a small scale and comparable with the non-relativistic velocities of the incoming states.

Figure 4 .
Figure 4. Thermally averaged Sommerfeld factors as a function of the time variable z = m χ /T for different color configurations of DVLF pairs.The relative mass splitting is ∆m/m χ = 10 −2 .

Figure 5 .
Figure 5. Self-energy diagram in pNRQCD with an initial scattering state in a color-octet (sextet) and an intermediate color-singlet (color-antitriplet) bound state.The imaginary part of the oneloop diagram on the left is responsible for the bound-state formation process (ψ ψ) p [8] → B n [1] + g, whereas the imaginary part of the right diagrams account for (ψψ) p [6] → B n [3] + g.

2 5 2 9 TeV mχ=3. 1 TeVFigure 7 .
Figure 7. Dark matter energy density Ω DM h 2 as a function of the dark matter mass m χ .The model parameters and the Majorana/Dirac options are indicated at the top of each plot.Magenta vertical lines serve as orientation to display the largest dark mater mass compatible with the observed energy density as obtained with a free cross section and with near-threshold effects (both Sommerfeld and bound states) respectively.

Ω DM h 2 5 Ω DM h 2 Figure 8 .
Figure 8. Dark matter energy density Ω DM h 2 as a funtion of the leptoquark mass m φ .The other model parameters indicated at the top of each plot, as well as the nature of the DM fermion (Majorana left panel, Dirac right panel).

Figure 9 .
Figure 9. Curves in the parameter space (m ψ , m φ ) that reproduces the observed DM energy density for y = 0.1, 0.5, 1.0, and for the Majorana and Dirac DM options respectively in the upper and lower panels.Two benchmark values for relative mass splitting are chosen, ∆m/m χ = 0.05 and ∆m/m χ = 0.01; the leptoquark-to-SM couplings are y L = y R = 1.The gray shaded area are the LHC limits on the masses of the DVLF and leptoquark of our model (for simplicity, here we quoted the most stringent bounds, which can be relaxed depending on the details, as explained later in the text).

Figure 10 .
Figure10.Left panel: The Dirac DM annihilation cross-sections into leptoquarks, χχ → φφ † with respect to the DM mass, m χ for the leptoquark mass, m φ = 1.5 TeV and two mass-splittings between DM and DVLQ, ∆m/m χ = 0.001 and 0.01, respectively, where all satisfy the correct relic abundance.Right panel: Photon spectra generated by different two-body and four-body final SM states, calculated using PPPC4DMID[84].

Figure 11 .
Figure 11.Region inside the red lines (1σ) corresponds to the required new-physics contributions to the (g − 2) µ via top-quark mass flip in the loop in the example model with φ = (3, 1, 1/3).For simplicity, here we quoted the most stringent collider bound (gray-shaded area), which can be relaxed depending on the model/analysis details, see discussion in the text.The purple band gives the observed DM energy density for m χ = 2.5 TeV and ∆m/m χ = 10 −2 .

96 Figure 12 .
Figure12.Region inside the red lines corresponds to the required new-physics contributions to the (g − 2) µ via charm-quark mass flip in the loop in the example model with φ = (3, 1, 1/3).The green-shaded area corresponds to the experimental result from NA62[93], whereas the gray-shaded area to non-resonant dilepton searches at the LHC[97,98].The purple band gives the observed DM energy density.