Leptoquarks and Real Singlets: A Richer Scalar Sector Behind the Origin of Dark Matter

We investigate scenarios with $\mathcal{O}(1\mathrm{\,TeV})$ scalar leptoquarks that act as portals between the Standard Model and Dark Matter. We assume that Dark Matter is a scalar singlet $S$ which couples to a scalar leptoquark $\Delta$ and the Higgs boson via the terms in the scalar potential. In addition, the leptoquark is endowed with Yukawa couplings to quarks and leptons that may address the anomalies in $B$ meson decays. We consider the $SS$ annihilation cross sections to estimate the Dark Matter relic abundance and explore the interplay between astrophysical, collider and flavour physics bounds on such models. In the heavy Dark Matter window, $m_S>m_\Delta$, the leptoquark portal becomes the dominant mechanism to explain the Dark Matter abundance. We find that the leptoquark Yukawa couplings, relevant for quark and lepton flavour physics, are decoupled from the dark matter phenomenology. By focussing on a scenario with a single leptoquark state, we find that relic density can only be explained when both $\Delta$ and $S$ masses are lighter than $\mathcal{O}(10\mathrm{\,TeV})$.


I. INTRODUCTION
Leptoquarks (LQs) are theoretically motivated hypothetical bosonic degrees of freedom that couple at tree-level to quark-lepton pairs [1,2]. They naturally appear in theories unifying quarks and leptons [3,4] and composite Higgs models [5,6], and they provide a viable mechanism to explain neutrino masses [7,8]. Their interactions with fermions also make them good candidates to address phenomenological problems in quark-lepton transitions, such as the discrepancies in B-meson decays observed at LHCb and the B-factories [9], or to address discrepancies in chirality-suppressed observables, such as the anomalous magnetic moment of leptons [10][11][12]. Whether these particles exist near the TeV scale remains a question to be answered by current and future experiments. In the meantime, it is natural to ask if LQs could also be related to other open questions in particle physics.
One of the most striking motivations for physics beyond the Standard Model (SM) is the evidence for Dark Matter (DM). Astrophysical and cosmological observations have accumulated indisputable * Electronic address:francesco.deramo@pd.infn.it † Electronic address:nejc.kosnik@ijs.si ‡ Electronic address:federico.pobbe@studenti.unipd.it § Electronic address:aleks.smolkovic@ijs.si ¶ Electronic address:olcyr.sumensari@physik.uzh.ch evidence at vastly different length scales [13]. These observations infer the presence of DM through its gravitational effects, but they tell us nothing about its microscopic nature. There are constraints from observations such as precise determination of its relic density [14], strong bounds on its electromagnetic interactions [15,16], and the fact that it must be stable on cosmological timescales [17,18]. Nevertheless, these constraints still leave room for a plethora of particle candidates [19].
An appealing scenario is the one where DM particles interact with the visible world enough to achieve thermal equilibrium in the early universe. As the universe expands and cools down, interactions between DM particles become less frequent. Thermal equilibrium is lost eventually, and such a departure from equilibrium is the physical process that sets the relic density. Remarkably, the resulting value depends only on dark sector masses and couplings that we can test in our experiments [20].
In this paper, we explore scalar LQs as mediators to the dark sector and whether this type of portal could explain the observed DM properties. We focus on scenarios where the DM particle is a scalar singlet that interacts with one of the scalar LQs via quartic interactions in the scalar potential. This scenario was the focus of Ref. [21]. We perform a comprehensive analysis that extends previous studies in several ways. We will consider the impact of all scalar LQ representations and provide the most general expressions for DM phe-nomenology. Moreover, we will study the interplay between the DM relic abundance with theoretical considerations of stability and perturbativity of the scalar potential, with direct and indirect DM detection constraints, as well as with flavour physics bounds. In particular, we will show that the LQ couplings needed to explain the observed DM relic abundance have little impact in flavour physics phenomenology. Furthermore, we will argue that the LQ and DM masses in the viable models must satisfy m ∆ < m S < O(10 TeV) in order to comply with the constraints from the stability and perturbativity of the scalar potential.
Several works have recently explored the possible connection between LQs and the observed DM abundance. Besides the possibility of scalar singlet DM considered in this work, the DM particle can be a fermionic singlet with gauge-invariant interactions with specific scalar and vector LQ representations [22,23]. DM can also belong to higher dimensional electroweak multiplets which can couple to the Pati-Salam vector LQ U 1 = (3, 1, 2/3) [24], which was proposed as a viable candidate to explain the so-called B-physics anomalies [25][26][27]. The connection between DM and the B-anomalies has also been explored in the context of composite vector LQ scenarios in Ref. [28].
Our setup and conventions are introduced in Sec. II. We list the phenomenological constraints on our framework in Sec. III. In particular, we compute the DM relic density, we impose experimental bounds from DM direct and indirect searches, and we consider constraints from collider physics. We investigate a concrete scenario with DM mass above the weak scale in Sec. IV. We discuss the possible connection with flavour anomalies in Sec. V, and we conclude in Sec. VI. Technical details are deferred to appendices.

II. A RICHER SCALAR SECTOR: LQS AND DM
In this Section we present our framework. We write down Yukawa interactions allowed by gauge invariance for different LQ representations, and scalar potential interactions for the DM candidate. We adopt the notation of Refs. [1,2] and specify LQ states by their SM quantum numbers (SU (3) c , SU (2) L , Y ).
We define the covariant derivative where Y is the hypercharge, and T A and T k are the relevant SU (3) c and SU (2) L generators, respectively. After the electroweak symmetry breaking, we can write where we idenfity the electric charge generator Q = Y + T 3 , the weak isospin raising/lowering matrices T ± = T 1 ± i T 2 , and the Weinberg angle θ W (we write for shortness s W ≡ sin θ W and c W ≡ cos θ W ).

II.1. Yukawa interactions between LQs and SM fermions
Scalar LQs have Yukawa couplings with the SM quarks and leptons. In order to describe the general features of LQ portals, we provide the general and convenient parameterization L Yuk = a,i,j q i y R ij P R + y L ij P L l j ∆ (a) + h.c. (3) These interactions are understood to be written for fermion mass eigenstates, and the Yukawa couplings y L and y R are matrices in flavour space. The chiral projectors P L,R isolate the corresponding chiralities of the fermion fields they act on. In our notation, l i stands either for a charged lepton or a neutrino. On the contrary, q i can be either a quark or its charge-conjugate. Finally, the scalar field ∆ (a) stands for a specific component (a) of the LQ multiplet ∆. We introduce concrete realizations for the interactions in Eq. (3). Following Ref. [1], we introduce the fermion number F defined as F = 3B + L, with B and L the baryon and lepton numbers, respectively. We assume there are no fermionic SM singlets.
Two scalar LQs can couple to fermion currents not carrying a net fermion number Left-handed fermion doublets in the above expressions are defined as where V and U denote respectively the Cabibbo-Kobayashi-Maskawa (CKM) and Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrices. Since neutrino masses are irrelevant for the phenomenology we will discuss, we are free to set U = 1. The Yukawa interactions in Eqs. (4) and (5) conserve the fermion number F as LQs have F = 0 in these cases. Baryon and lepton numbers are conserved upon assigning them to the LQ field as well.
Alternatively, one can consider LQ interactions with fermion currents carrying a net |F | = 2 fermion number where Ψ C denotes a charge-conjugated fermion. These LQ states are potentially dangerous because they could have diquark couplings which would trigger baryon and lepton number violation [2,29].
Here and after, we assume that these couplings are forbidden by a suitable symmetry that guarantees proton stability. Finally, we provide a prescription to connect the generic couplings y L ij and y R ij defined by Eq. (3) onto ones of the specific models listed in Eq. (4)- (8). This matching is given in Tabs. I and II for charged leptons and neutrinos, respectively.

II.2. Higgs and LQ portals to Dark Matter
We extend the framework introduced above by adding the DM candidate, a real scalar S not carrying any SM gauge quantum numbers. Stability of S is ensured by a Z 2 symmetry under which S → −S while other fields remain unchanged. Interactions between dark and visible sectors, at the renormalizable level, proceed via scalar potential couplings. Here, we provide the phenomenological expressions for these interactions which can be viewed as a generalization of the Higgs portal model [30][31][32].
The general Lagrangian for the scalar sector takes the form where H denotes the SM Higgs doublet field and we consider a single LQ multiplet ∆. Besides the kinetic terms, with covariant derivatives as given in Eq.
(1), we have a scalar potential that we parameterize as follows The first contribution is the same one as in the SM involving just the Higgs doublet Here, we have µ 2 > 0 in order to ensure the spontaneous breaking of the electroweak symmetry. The BSM scalar potential contains masses and interactions of the new degrees of freedom where the free real parameters m 1,2 and λ i are always allowed by the symmetries of the theory. For the case of LQ weak doublets (i.e., for ∆ = R 2 or R 2 ), the additional term is allowed by electroweak gauge invariance. Once the SM Higgs gets a vacuum expectation value (vev), this interaction induces a mass splitting among the SU (2) doublet LQ states. This operator, also constrained by T -parameter [33], is not

Neutrinos
Coupling q l R2 R2 S1 S3 S1  [34]. These operators can be relevant to generate neutrino masses [35], or to generate dipole operators to accommodate the muon g − 2 anomaly [36]. Nonetheless, they do not play any relevant role in DM phenomenology. For these reasons, we neglect these further interactions in our study.
Not all couplings in the scalar potential in Eq. (12) have the same impact on our analysis. The quartic couplings λ 1,2 , which correspond to self-interactions of S and ∆, do not play a considerable role in phenomenological analysis apart from their impact on the stability of the potential. Focusing on interactions that do impact DM phenomenology, we can express the scalar potential after electroweak symmetry breaking as follows Here, h is the SM Higgs boson and v = ( √ 2G F ) −1/2 is the Higgs vev. DM phenomenology will change accordingly to the values assumed for the quartic couplings λ 3,4,5 , and the masses m S and m ∆ .

II.3. Stability constraints
Here we discuss the theoretical constraints on the quartic couplings λ i arising from the stability of the scalar potential in Eq. (10). The quartic terms in Eqs. (11) and (12) can be written as a quadratic form The above expression has to be bounded from below for large field values in any arbitrary direction. If we consider going to infinity along the direction of a single field, |H| 2 → ∞, S 2 → ∞, or |∆| 2 → ∞, this results in necessary conditions λ, λ 1 , λ 2 > 0, respectively. In order to find the sufficient conditions, we diagonalize the n × n matrix appearing in Eq. (15) and require all its eigenvalues x 1,2,...,n to be positive. This quadratic form is diagonalizable and thus its characteristic polynomial p(x) has n roots x 1,...,n and n − 1 stationary points e 1,...,n−1 which are ordered as If some of the roots are multiple, then the above inequalities become equalities across the range spanned by multiple roots (e.g., for a double root We show in Fig. 1 the characteristic polynomial p(x) for the case when all roots are positive. As we can see, the requirement x 1 ≥ 0 is fulfilled precisely when e 1 ≥ 0 and p(0) ≥ 0 since p(x) > 0 for large negative x. If we repeat the above argument for e 1 ≥ 0, it follows that p (0) ≤ 0 and that the first stationary point of p (x) should be non-negative. Iterating this argument to yet higher derivatives yields the conditions 1 Thus sufficient stability conditions can be expressed in terms of λ i . We find it convenient to define the parameters 1 This is a special case of the so-called Descartes' rule of signs.
In addition to the already mentioned positivity of λ and λ 1,2 , we derive the following sufficient stability conditions |x|, |y|, |z| < 1 , The inequalities (19) define a tetrahedron-like region in the {x, y, z} space. Are the above conditions also necessary? It turns out that since the potential matrix is multiplied by strictly positive vectors from the left and right, the above condition can be somewhat relaxed. We require that matrix V is co-positive definite, which by definition means that x T V x > 0 for any vector x from the first octant. Following Ref. [37], in particular Eqs. (5) and (6) of that paper for particular case of 3 × 3 matrices, we find slightly looser conditions compared to (19): Together with λ, λ 1 , λ 2 > 0, these are necessary and sufficient stability conditions. The most interesting among these inequalities are the ones for λ 3,4,5 since they are related to the interactions between the scalar particles in our model. In particular, we can derive the following necessary conditions on λ 3,4,5 , where the leftmost bounds are realised when both λ 1 and λ 2 are set to the perturbativity limit λ 1,2 = √ 4π and inserted to the stability lower bounds. The rightmost upper bounds stem from the perturbativity requirement, whereas stability considerations do not yield an upper bound on λ 3,4,5 .

III. HIGGS AND LEPTOQUARK PORTALS PHENOMENOLOGY
In this section, we study the phenomenology of our scalar leptoquark portal framework. DM phenomenology is driven by the scalar potential interactions in Eq. (14), and the scalar singlet S interacts with the visible world via the Higgs portal (λ 4 coupling) as well as through the LQ portal (λ 3 coupling). In this section, we compute the DM relic density, we study constraints from direct and indirect searches and we account for collider bounds.

III.1. Relic density
The DM number density n S evolves according to the Boltzmann equation [38,39] Here, the Hubble parameter H = 1 a da dt accounts for dilution due to the expansion and it is defined in terms of the cosmological scale factor a and its derivative with respect to the cosmic time t. The change in the number of S particles due to collisions is accounted for by the term on right-hand side, proportional to the thermally averaged annihilation cross section times the Møller velocity σv . The inverse process, where degrees of freedom from the thermal bath annihilate to produce a pair of DM particles, is proportional to the square of the DM equilibrium number density n eq S .
When DM particles are in equilibrium in the early universe, their energies and momenta are distributed in phase space accordingly. Thus the kinematics of the initial state for each DM annihilation is not fixed, and we need to account for all possible options. This is the reason why we have a thermal average in Eq. (22). There is a general expression connecting the annihilation cross section as a function of the Mandelstam variable s, which is nothing but the (square of the) center of mass energy, and the thermally averaged cross section. Here, T is the temperature of the primordial thermal bath and K 1,2 (x) are modified Bessel functions of the second kind. This expression is valid for a Maxwell-Boltzmann statistics of initial state DM particles, and we can trust it since quantum degeneracy effects give very small corrections in the early universe. In order to determine the relic density, we need to compute expressions for all the DM annihilations allowed by kinematics. The phenomenologically relevant DM annihilation channels to visible final states are: (i) SS → hh, (ii) SS → ∆∆, (iii) SS → V V (gauge bosons) and (iv) SS → f f (fermions). We account for the leading contributions for each channel, which can appear either at tree or one-loop level and depend on the underlying parameters of the scalar potential as well as leptoquark flavour parameters. These expressions should be convoluted with the thermal distributions at finite temperature T as prescribed by Eq. (23). We provide full derivations and complete expressions for the annihilation cross section in App. A and we evaluate the DM relic density by following standard techniques [38].
The leading annihilation channel, among the several ones available, depends on what parameter space region we focus on. For illustration, we consider two different benchmarks where we fix the quartic couplings λ 3,4,5 . We neglect for this discussion the LQ Yukawa couplings to SM fermions; we comment their impact on DM phenomenology in Sec. IV.
Higgs portal: λ 3 = 0, λ 4,5 = 1: The leading annihilation channels are driven by the λ 4 quartic coupling with the SM Higgs field. The λ 5 quartic interaction between LQ and Higgs field plays only a marginal role since it only contributes to the subdominant channel SS → h → ∆∆ * proportionally to λ 4 λ 5 .
We show the thermally-averaged annihilation cross sections as functions of the DM mass m S for this scenario in Fig. 2. For the purpose of illustration, the LQ state ∆ is in the R 2 = (3, 2, 7/6) representation of the SM gauge group. We provide three snapshots at three different temperatures T /m S ∈ {1, 1/20, 1/50}; this is the relevant temperature range for thermal freeze-out of cold relics. Annihilations via s-channel Higgs exchange are by far dominant. The dominant final states can be either ff or V 1 V 2 (with V 1,2 a SM electroweak gauge boson), depending on the mass of S, whereas annihilations to LQs are sub-dominant since LQs do not couple directly to DM in this scenario. For this reason, LQ interactions only contribute indirectly to the relic density via the modification of the hV 1 V 2 couplings, which are mostly relevant for m S values around the electroweak scale. The thermally averaged cross sections have a peak at the Higgs pole, m S = m h /2, which becomes wider but still dominant as the temperature T increases. We comment about the (in)viability of this scenario to ac-   LQ portal: λ 3 = 1, λ 4,5 = 0: In this case, the situation is inverted since DM does not couple at tree-level to the Higgs boson whereas it has O(1) couplings with LQs. Once we set λ 4 = 0, this scenario is actually completely blind to λ 5 since it turns out that only a combination of λ 4 λ 5 enters the annihilation cross sections. From Fig. 3, we see that the dominant annihilation channel is to LQ final states, which become kinematically accessible above the threshold m S > m ∆ , via the λ 3 quartic interaction. For smaller values of m S , the dominant modes are SS → V 1 V 2 induced by LQ loops. Cross sections in Fig. 3 are much smaller than the ones appearing in Fig. 2 since they are loop suppressed, except for the ∆∆ * channel, and since they are not resonantly enhanced as in the case of the Higgs portal regime. Here, the relevant parameter space will appear at large DM masses where direct detection constraints are less effective, as we discuss in Sec. IV.
The above discussion was based on the R 2 representation for the LQ state. Considering different representations would only affect the annihilation channels into the gauge bosons V 1 V 2 , since they depend on the SU (2) L and U (1) Y quantum numbers of the LQs running in the loops, see e.g. Eq. (B7) in App. B, which amount to mild modifications of the annihilation cross section. On the other hand, the ∆∆ * channel is enhanced for LQs with larger multiplicity 2T + 1, but without changing the general features described above.  [40] in the plane mS vs. |λ3| (left panel) and mS vs. |λ4| (right panel). We fix m∆ = 1.5 TeV and the other couplings are set to zero. Constraints on λ3 depend on the specific SU (2)L representation (singlet, doublet or triplet), while the ones on λ4 are, to first approximation, independent of LQ interactions. The coupling λ5 cannot affect direct detection constraints (see text).

III.2. Direct detection
DM effective interactions with nucleons N mediate elastic scattering that can be searched for by direct detection experiments. The effective operator relevant to our analysis is where C N is a low-energy Wilson coefficient. From the microscopic point of view, these interactions arise from DM couplings to quarks and gluons. The corresponding Lagrangian, defined at the scale µ = O(1 GeV) above confinement, is [41] L eff QCD = q=u,d,s Only effective couplings to light quarks, C q , appear in the Lagrangian since heavy quarks (c, b, t) have been integrated out. Their virtual effects provide further contributions, besides the ones due to new physics, to the effective couplings to gluons C g . We can match these effective coefficients to quarks and gluons onto the nucleon Effective Field Theory [42] The quantities f N T q and f N T g are set by nucleonic matrix elements. We use the following values: .032 and f p Ts = 0.020 for a proton and f n Tu = 0.017, f n T d = 0.041 and f n Ts = 0.020 for neutrons, with f p,n T G = 1 − q=u,d,s f p,n T q [43]. This operator induces a spin-independent scattering between the DM particle and the target nucleus. We report the cross section for this scattering normalized per nucleon where Z and A are the atomic and mass numbers of the nucleus, and µ p = m p m S /(m p + m S ) is the reduced mass for the DM-proton system. The effective low-energy couplings for our LQ scenario can be read from Eq. (B13) where q ∈ {u, d, s} and n H = 3 is the number of heavy quarks (H). Contributions of order 1/m H in the heavy-quark expansion have been neglected. These expressions extend results already existing in the literature [21], and the new contributions we compute give order one corrections to the scattering cross section. In particular, the contribution proportional to λ 5 gives a small correction due to the additional suppression factor v 2 /m 2 ∆ .
The current most stringent bounds on spinindependent DM-nucleon cross section come from the XENON1T experiment [40]. In Fig. 4, we show the constraints on |λ 3 | (left panel) and |λ 4 | (right panel) as a function of m S , for a benchmark value m ∆ = 1.5 TeV. The most constrained coupling turns out to be λ 4 , independently of the LQ mass, and we find that the dominant contribution comes from the Higgs gluonic penguins.

III.3. Indirect detection
DM annihilations produce final state photons that are searched for with gamma-ray telescopes. The galactic center (GC) is a natural target and it comes with pros and cons. On one hand, it is not excessively far away and we do not lose too much flux. On the other hand, it is a very active region and astrophysical backgrounds are significant and not completely understood. Dwarf spheroidal satellite galaxies (dSphs) of the Milky Way are immune from this issue since they are DM-dominated and with considerably fewer backgrounds. In our study, we consider bounds from both of these sources. The strongest constraints for DM masses around or below the weak scale come from the Fermi Large Area Telescope (Fermi-LAT) [44][45][46]. Different instruments provide meaningful bounds on TeV scale DM candidates: the High Energy Stereoscopic System (HESS) observatory [47][48][49], the Very Energetic Radiation Imaging Telescope Array System (VERITAS) [50], and the Major Atmospheric Gamma Imaging Cherenkov Telescopes (MAGIC) [51]. For even larger masses, bounds from the High Altitude Water Cherenkov (HAWC) Observatory [52] are the most severe ones. Besides imposing current experimental constraints, we also exploit the discovery reach of the future Cherenkov Telescope Array (CTA) [53].
Experimental collaborations provide bounds on the annihilation cross section as a function of the DM mass for fixed SM final states. There are several annihilation channels available in our framework, and we identify the dominants ones for the two benchmarks introduced at the beginning of this Section: the Higgs portal and the LQ portal. Higgs and weak gauge boson final states are the dominant channels for the former case. Performing an analysis by computing the gamma-ray spectrum from all different final states is beyond the scope of our work, and we provide conservative bounds by imposing the most severe constraints among the ones for the final state particles; our indirect detection bounds are thus rather conservative for this bench-mark. For the LQ portal case, we have two LQs in the final state that decay subsequently to quarks and leptons; the DM annihilate to four final state SM particles [54]. We compute the expected photon spectrum from annihilations to LQs and we compare with the experimental sensitivities of current and future experiments. As we discuss in the next Section, the region of our interest is the one where the DM mass is above the weak scale. We can establish the dominant annihilation channels for our benchmarks by looking at the right panel of Figs. 2 and 3; DM particles are non-relativistic today and therefore annihilations are correctly described by thermal averages at low temperatures.

III.4. Collider constraints
We conclude this phenomenological section with a list of collider constraints. Conventional mono-X searches for DM do not provide the most stringent bounds from collider physics for our framework. The scalar potential couplings of the Higgs doublet H to S and ∆ leave an imprint on Higgs properties measured at the LHC. Thus we consider Higgs physics as well as direct searches for LQs.

III.4.1. Invisible Higgs Decay
For a DM particle coupled to the Higgs boson which is light, m S < m h /2, an important constraint comes from the upper bound on the branching fraction of Higgs decays into invisible final states [55][56][57][58][59][60]. Within our framework, we have the requirement where Γ SM h = 4.1 MeV [61]. The invisible partial decay width results in We impose the experimental bound Br(h → inv.) < 0.24 [62] and we find |λ 4 | 0.02 for m S values below 50 GeV. The constraint is significantly relaxed when m S is just below m h /2, according to phase space suppression in Eq. (31). This constraint excludes the viability of the low-mass DM since the annihilation cross section of SS, driven by the Higgs portal, λ 4 , is too small which leads to overabundance of S.

III.4.2. Higgs Decay to Leptons
The measurements of h → µ + µ − and h → τ + τ − provide important constraints on the µt and τ t leptoquark Yukawas via the modification of the Higgs Yukawa in Eq. (A12). Recently, the CMS collaboration 2 provided the best-fit values of the signal strength parameters as [64] For the τ τ channel, the CMS collaboration also provided [65] Using these constraints and assuming that the modification on the Higgs production due to LQs is negligible, we obtain the following 1σ constraints from Eq. (A12) of App. A on products of the LQ Yukawa couplings,

III.4.3. Higgs production at the LHC
The production cross section of the Higgs boson at the LHC via gg → h and its decay h → γγ are affected by the virtual LQ loops, proportional to the λ 5 coupling. In the approximation where these are the dominant effects of leptoquark loops, we can use the combined result of ATLAS and CMS (Fig. 17 of Ref. [66]), which constrains the relative SM coupling modifiers κ g and κ γ . Explicit expressions for κ g,γ are where x i = m 2 h /(4m 2 i ), as in Eqs. (B17) and (B18) in App. B. We find that the current experimental precision on κ g,γ − 1 0.2 ≈ λ 5 v 2 /(4m 2 ∆ ) results in rather weak constraint λ 5 /(4m 2 ∆ ) 3/TeV 2 .

III.4.4. LHC direct searches
LQs can be produced in pairs at the LHC via gluon fusion and decay into quark-lepton pairs. Searches have been performed at ATLAS and CMS for several possibilities of final states, cf. e.g. [67][68][69] for a compilation of the latest LHC bounds. For instance, for LQs decaying into tµ final states, we find that the most stringent upper limits on LQ masses are 1420 (950) GeV for B(∆ → tµ) = 1 (0.5) [70]. Similar or weaker limits are obtained for other possible final states. In what follows, we will conservatively take m ∆ ≥ 1.5 TeV, allowing us to safely avoid existing LHC direct-search limits.

IV. A VIABLE SCENARIO WITH HEAVY DARK MATTER
We assess the viability of possible DM scenarios by combining the phenomenological constraints presented in the previous Section. For DM lighter than the weak scale, severely constrained by XENON1T as illustrated in Fig. 4, annihilations to SM fermions via the Higgs boson exchanged in the s-channel are dominant. This case is potentially interesting because LQ Yukawas, motivated by discrepancies observed in B-meson decays (see e.g. Ref. [68,71] and references therein) and in the anomalous magnetic moment of the muon [36,72,73], could alter the naïve expectation for the annihilation rates. More precisely, these LQ flavour effects would come into play via modification of the Higgs Yukawas contributing to SS → h * → via a chiral enhancement of the amplitudes (∝ m t /m ) which lifts up the cross section significantly compared to the standard Higgs portal.
In order to show that the low-mass DM regime is not viable, we take the LQ R 2 = (3, 2, 7/6) as benchmark since the chiral enhancement of the leptonic amplitude is possible for this model. If we set the LQ Yukawa couplings to zero, the quartic λ 4 controls annihilations to SM fermions in the low DM mass region, as depicted in Fig. 2. This is also true if we switch on the Yukawa couplings, with a chirally enhanced contribution to the amplitude for annihilation to lepton pairs + − proportional to y t L y t R (see Eq. (A12) of App. A). We look for the  Figure 5. The values of the Higgs portal coupling λ4 needed to explain the DM abundance are plotted in blue as a function of the DM mass mS. Results for singlet (S1, S1), doublet (R2, R2) and triplet (S3) LQ states are shown, with m∆ = 1.5 TeV. These results are largely independent on the hypercharge of the LQ state, depending only on the SU (2)L representation. Excluded regions due to Xenon1T and perturbativity are depicted by the gray-shaded areas. The blue-shaded regions are excluded by indirect constraints by HESS [47] and FermiLAT [45]. Note that in the right-most plot there are forbidden intervals of mS, as there the annihilation cross-section is too large due to the large value of λ3.
values of λ 4 necessary to explain the observed relic density for different values of |y t L y t R | within the perturbative regime, and we compare them with bounds on λ 4 stemming from the invisible Higgs width in Eq. (30) and the XENON1T limits depicted in Fig. 4. In the low DM mass region, there is no value of λ 4 compatible with relic density and not excluded by the constraints mentioned above. The couplings y L t y R t provide a non-neglibible decrease in the required value of λ 4 but their effect cannot bring λ 4 to the O(10 −2 ) level which is needed.
In order to be compatible with the severe experimental constraints at low DM mass, we consider values of m S above the electroweak scale. LQ Yukawas are irrelevant in this DM mass region since they do not impact the main annihilation channels illustrated in Fig. 2 and 3, and annihilations are controlled by the quartic couplings λ 3,4,5 . In what follows, we explore in detail the viability of such a heavy DM scenario.
We now consider the singlet (S 1 , S 1 ), doublet (R 2 , R 2 ) and triplet (S 3 ) LQ states ∆ with mass m ∆ = 1.5 TeV, in agreement with the LHC direct limits from Sec. III.4. The leading impact of the LQ hypercharge (Y ) is to modify the subleading annihilation channel SS → ZZ, where its contribution is suppressed by s 4 W (cf. Eq. (B10)). For this reason, we observe no noticable effect of Y in Figs. 5 and 6. Furthermore, we neglect the LQ Yukawa couplings since they have negligible impact on DM phenomenology in this mass region. We show in To further explore the regime where SS → ∆∆ * acts as the main annihilation channel, we plot in Fig. 6 the allowed DM parameters in the (m S , λ 3 ) plane by keeping λ 5 = 0 and fixing this time λ 4 ∈ {0, 0.05, 0.5}. For the pure LQ portal (λ 4 = 0) shown in the left panel, we see that a coupling λ 3 0.5 is needed to reproduce the DM relic abundance. For larger values of λ 4 , both SS → ∆∆ * and SS → V V annihilation channels become relevant, and smaller values of λ 3 are needed as a consequence. Direct detection constraints are mostly sensitive to λ 4 and practically insensitive to λ 3 . There are two potential constraints on the LQ portal coupling λ 3 : indirect detection and perturbativity. The former is sub-dominant, as we quantify in the next paragraph, and the only meaningful con- Figure 6. The values of the LQ-portal coupling λ3 needed to explain the DM abundance are plotted in blue as a function of the DM mass mS. Results for singlet (S1, S1), doublet (R2, R2) and triplet (S3) LQ states are shown for m∆ = 1.5 TeV. DM constraints turn out to be largely independent on the hypercharge of the LQ state, depending only on the SU (2)L representation. Excluded regions due to Xenon1T and perturbativity are depicted by the gray-shaded areas. Note that below the ∆∆ * threshold the value of λ3 becomes very large in the first two plots, whereas in the last plot (λ4 = 0.5) there is no solution below a certain mass mS close to the threshold since the annihilation cross-section is too large.  Figure 7. Flux of gamma rays from DM annihilations to LQs and subsequent LQ decays to SM particles (dashed orange lines). We consider both LQ decays to third generation fermions (tτ ) as well as lighter fermions (qe). The annihilation cross section reproduces the relic abundance via thermal freeze-out, and we consider a dwarf galaxy with the J-factor provided in the legend. We report for comparison also the expected spectrum for a WIMP candidate annihilating to bb pair for different masses (dotted black lines). The sensitivity lines for current (solid) and future (dashed) experiments are taken from reference [74]. straint on λ 3 arises from the perturbative bounds.
We compute the gamma-ray spectrum produced by the cascade annihilations, first with DM anni-hilating to LQ final state SS → ∆∆ * and then subsequent LQ decays, and we show our results in Fig. 7. We fix the annihilation cross section to the Figure 8. Constraints on the LQ-portal λ3 and the Higgs-portal λ4 couplings imposed by the relic abundance (reproduced exactly along the blue line), direct detection due to Xenon 1T (exclusions with a dashed border) and the potential stability constraints (exclusion with a dash-dotted border). Results for singlet (S1, S1), doublet (R2, R2) and triplet (S3) LQ states are shown, with two representative values of mS above m∆ = 1.5 TeV. DM constraints turn out to be largely independent on the hypercharge of the LQ state, depending only on the SU (2)L representation.
thermal relic value, and we consider a dwarf galaxy with a J-factor equal to J = 10 19 GeV 2 cm −5 . The dashed orange lines provide the photon flux for the representative spectrum (m S , m ∆ ) = (5, 1.5) TeV and for two possible LQ decays: third generation fermions (∆ → tτ ) and lighter fermions (∆ → qe). The larger top mass is responsible for a slight enhancement in the UV tail of the spectrum. We report for comparison also WIMP spectra for annihilations to bb final state for different DM masses. The only region where the LQ portal spectrum differs slightly from the one induced by annihilations of a WIMP with the same mass is in the UV tail. We show in the same figure the sensitivity curves for different instruments, as collected in Ref. [74]. The Fermi bounds for DM mass around the TeV scale in the LQ portal benchmark are for all purposes the same as the one for a WIMP with the same mass; thermal relics are not excluded. For instruments reaching maximum sensitivity around the TeV scale, the slight enhancement in the UV tail could slightly affect the bounds for WIMPs but they should not spoil the picture completely; these instruments are still quite far from the thermal relic line and therefore thermal relics are also not excluded by them. We conclude that current indirect detection constraints are never the most stringent ones for the LQ portal case. We also compare the sensitivity of the future CTA telescope with the signal predicted in our model.
We visualize the interplay among the different constraints in Fig. 8 where we fix m S to two benchmark values in the TeV range, and we vary the LQ and Higgs portal-couplings λ 3 and λ 4 . We notice how an improvement of the direct-dection constraints by a factor of ≈ 4 is needed to start probing λ 3 in this range of DM masses. Moreover, stability constraints (Eq. (20)) turn out to be meaningful for negative values of λ 4 .
Finally, we explore the extent to which our conclusions remain valid when the LQ mass is increased above our benchmark value m ∆ = 1.5 TeV. To this purpose, we focus on the pure LQ portal, with λ 4 set to zero. The values of λ 3 needed to reproduce the relic density are then shown in Fig. 9 in the plane m ∆ vs. m S for singlet, doublet and triplet LQ states. The boundary of the allowed region is determined by the perturbativity constraints, which imply the following upper limit on both the LQ and DM masses The perturbativity contours for the doublet and triplet LQ states are broader than the correspond- ing singlet contour. This is a consequence of higher LQ multiplicity, which increases the annihilation cross section and slightly relaxes the upper limit Eq. (38). However, the main phenomenoligcal features remain very similar for all scenarios.

V. IMPLICATIONS FOR FLAVOUR ANOMALIES
In this Section, we discuss the implications of the viable DM scenarios outlined in Sec. IV to the discrepancies observed in B-meson decays and the muon g − 2. Both types of discrepancies suggest that LQs could exist at the TeV scale with O(1) Yukawa couplings to the SM fermions. Therefore, it is a natural question if these anomalies could be accommodated in a scenario that would also explain the DM abundance via the mechanism discussed in Sec. IV.
To answer the question raised above, we will remind the reader which scalar LQ states can successfully accommodate each of these anomalies and we will derive the ranges of masses m ∆ that are compatible with them [75,76]. These values will be compared to the upper limits on m ∆ obtained in Fig. 9 in order to provide a viable DM scenario. As already anticipated in Sec. IV, the specific values of the LQ Yukawas that are fixed at low-energies have little impact on the DM abundance, in such a way that the only common parameter for flavor and DM phenomenology is indeed the LQ mass, which is explored in the following.
integrated in the di-lepton invariant-mass bin q 2 ∈ (1.1, 6) GeV 2 , turns out to be 3.1 σ below the precise SM prediction, R SM K = 1.00(1) [78,79]. Deviations from the SM predictions have also been observed in similar LFU tests with B → K * decays [80]. The combination of these results with the current experimental average of the B s → µµ branching fraction [81, 82] which is also slightly below the clean SM prediction B(B s → µµ) SM = 3.66(14) × 10 −9 [83], amounts to a combined deviation of 4.6 σ from the SM predictions [69] (see also Ref. [84][85][86]). The preferred scenario to explain these deviations requires a purely left-handed operator, with the effective coefficient C µµ L in the following range [69], Among the scalar LQs listed in Sec. II, only the scalar triplet S 3 = (3, 3, 1/3) can induce a nonzero value of C µµ L at tree-level [69,87], where we remind the reader that the Yukawas y L S3 are defined in Eq. (7). By combining Eq. (42) and (43), it is straightforward to conclude that This constraint implies an upper bound on m S3 which is compatible with, but much less strict than the bound obtained in Fig. 9 from the requirement of reproducing the correct DM relic density via the mechanism discussed in this paper. The bound on m S3 from relic density becomes significantly stronger for small mass m S .
• R D ( * ) : Several discrepancies from the SM predictions have also been observed in LFU tests for the transition b → c ν, The current experimental averages of LHCb [88,89] and the B-factories [90][91][92][93][94] measurements are [95] R exp which turns out to be ≈ 1.3 σ and ≈ 3.5σ above the SM predictions, which are obtained by combining the latest lattice QCD results for the B → D ( * ) [96,97] form factors at nonzero recoil with the B → D ( * ) lν (l = e, µ) differential decay rates measured experimentally, see [97,98] and references therein. The observed deviations in R D ( * ) can also be interpreted by means of an effective field theory, where g a are effective couplings, defined at the renormalization scale µ = m b , and the relevant effective operators are as well as O V R and O S R which are obtained from the ones above by flipping the chirality of the quark fields. Differently from the previous case (see Eq. (41)), there is more than one viable effective scenario induced by LQs that can explain the anomalies in R D ( * ) . The simplest viable scenario is to consider once again a purely left-handed operator, where we have updated results from Ref. [69] by considering the latest lattice QCD results for the B → D * form factors [97]. Another option is to consider the effective scenarios g S = ±4 g T , at the matching scale µ = m ∆ ≈ 1 TeV, which can also perfectly describe current data [69]. These relations become g S L (m b ) ≈ −8.5g T (m b ) and g S L (m b ) ≈ +8.1g T (m b ) at µ = m b , respectively, after accounting for the RGE effects [99]. These scenarios can provide a viable explanation for the b → cτν anomalies for real and purely imaginary couplings, respectively [69] (see also Ref. [100,101]). There are only two scalar LQs that can predict the allowed values of effective couplings at lowenergies, while being consistent with various flavour and LHC constraints [69]. The first scenario that can be matched to the viable effective scenarios described above is S 1 = (3, 1, 1/3), which can predict nonzero values for both g V L and g S L = −4g T via the products of couplings y L S1 y L * S1 and y L S1 y R * S1 , respectively. The needed couplings to explain the anomalies can then be either y L S1 bτ V y L S1 * cτ or y L S1 bτ y R S1 * cτ which both require LQ massses in the O(TeV) range.
The second viable scenario is R 2 = (3, 2, 7/6), which can explain current data provided there is an imaginary phase in the LQ Yukawas [69,100], which is also in the O(TeV) range. Therefore, by inspecting Eq. (53)- (54) and (55), we find that the range of LQ masses needed to explain R D ( * ) turns out to be very similar to the one needed to explain the DM relic abundance, see Fig. 9. Note, also, that these masses and LQ couplings are fully consistent with high-p T constraints and with other flavor constraints, as recently analyzed in Ref. [69].
• (g − 2) µ : Lastly, we discuss the impact of scalar LQs to the 4.2 σ discrepancy observed between the experimental determinations of a µ = (g − 2) µ /2 [10,11] and the SM prediction from the Muon g − 2 theory initiative [12], This discrepancy is comparable in size to the SM electroweak contributions. Therefore, it is only possible to explain it through LQs contributions with O(TeV) masses if a chirality-enhancement mechanism takes place [102][103][104][105]. Such an enhancement can be induced by LQs that couple simultaneously µ L and µ R to a heavy fermion, which is typically the top quark, producing an enhancement ∝ m t /m µ [72]. 3 There are only two LQ states that satisfy this criteria, namely S 1 = (3, 1, 1/3) and R 2 = (3, 2, 7/6). The needed couplings for S 1 read y L S1 tµ y R S1 * tµ whereas for R 2 , where, for simplicitly, we have kept only the chirality-enhanced contributions and set the LQ masses to 10 TeV in the logarithms. In this case, we find that reproducing the observed DM relic density imposes a much stricter constain than the muon g − 2 which is sensitive to very large LQ masses. Note, also, that the constraints obtained in Eq. (57) and (58) are much stricter than the ones derived from other loop observables such as Z → µµ [107,108]. Finally, it is worth stressing that even though S 1 and R 2 are also the scalar LQs needed to explain R D ( * ) , the simultaneous explanation of both anomalies is tightly constrained by τ → µγ which would also be chirality-enhanced for the needed pattern of Yukawas [109].
To summarize this discussion, DM phenomenology and flavour physics are complementary probes of the LQ parameters. DM constraints are practically insensitive to the LQ Yukawas, but they can be used to fix the scalar potential parameters m S , m ∆ and λ 3,4,5 , whereas flavour-physics constraints, at tree-level, are only sensitive to the combinations of Yukawas |y ij |/m ∆ , where y ij denotes a generic LQ Yukawa with flavour indices i, j. For this reason, the only connection between flavour and DM arises from the LQ masses, which are bounded in both cases by perturbativity constraints, as shown e.g. in Eq. (38) for the quartic couplings and discussed above for each of these anomalies. By comparing these upper bounds with Eq. (38), we see that the limit on m ∆ derived from DM relic abundance is more constraining than the ones derived from most flavour anomalies, with the only exception of R D ( * ) , which is a tree-level process in the SM. Therefore, if any of these anomalies is confirmed in the future, one could easily check from Eq. (38) if the LQ scenario favored from flavour physics could be extended to accommodate DM via the mechanism discussed in this paper.

VI. CONCLUSION
The origin of DM is a long standing problem in physics of fundamental interactions. The recently observed discrepancies in the lepton flavour universality tests in B meson decays suggest the existence of LQs. In this paper, we have studied a possible connection between these two open problems. We took a real scalar singlet field S as a DM candidate, and introduced a scalar LQ ∆ with mass above 1 TeV and O(1) Yukawa couplings to the SM fermions which can resolve the flavour anomalies. The only renormalizable interactions of the DM field S with the visible fields are the Higgs and LQ portal-couplings present in the scalar potential.
We have studied the role of the LQ in DM phenomenology, focussing in particular on the portal coupling λ 3 |∆| 2 S 2 , mass m ∆ , and Yukawa couplings. The latter two parameters also play a central role in flavour anomalies. First, we have introduced the most general scalar potential for the scalar fields of our framework, and we have carefully studied the scalar potential stability conditions which yielded correlated upper bounds on the portal couplings. For the LQ Yukawa sector we stated explicitly all possible representations along with their couplings to the SM fermions. Assuming one scalar LQ is present, we have calculated the DM annihilation rates that depend on the scalar potential parameters and Yukawa couplings. Expressions for cross sections for general LQ representation are provided in App. A. We feed the Boltzmann equation for the S number density with these cross sections and compute the relic density. Up-to date constraints from direct and indirect searches for DM are also derived and they place important upper bounds on the scalar potential couplings. Collider constraints on S have been taken into account, including the very stringent upper bound on Higgs boson invisible width.
Two scalar potential couplings drive the DM phe-nomenology -the Higgs portal λ 4 |H| 2 S 2 and the LQ portal λ 3 |∆| 2 S 2 . It was found and shown in Fig. 5 that below the threshold for annihilation into LQs, i.e. m S < m ∆ , the loop effects of λ 3 are loopsuppressed and cannot compete with the Higgs portal coupling λ 4 . On the other hand, we could observe the impact of large LQ Yukawas boosting SS → annihilation cross sections ( = µ, τ ). This effect is potentially large at low DM mass where channel drives the total annihilation cross section, but the Higgs invisible width is too strong a constraint in this mass region and the relic density comes out too small, despite Yukawa enhancement. At higher m S the Yukawa couplings do not affect the S annihilation cross section significantly.
When m S > m ∆ the SS → ∆∆ * channel is open and the phenomenology depends on both Higgs and LQ portals. It turns out that large LQ portal λ 3 implies small Higgs portal λ 4 and vice versa. Although both scenarios are equally viable the Higgs portal regime is more prone to direct detection constraints as well as to the stability constraints as shown in Fig. 8. Further analysis of the pure LQ portal regime with λ 4 = 0 reveals a wide parameter space of λ 3 and m S > m ∆ 1.5 TeV. Most importantly, the requirements of stability of the scalar potential and perturbativity set limits both on m S and m ∆ to lie below O(10 TeV), as shown in Fig. 9.
In summary, flavour and DM aspects of the scalar singlet and scalar LQ model are to a large degree decoupled when m S is below the threshold for annihilation into LQs, whereas above the threshold we may enter the LQ portal regime where both S and ∆ masses are bounded from above due to the pertubativity constraint. The obtained mass bounds on the LQ are stricter than what would be inferred from most flavour anomalies.
(i) SS → hh : Feynman diagrams for annihilation into Higgs boson pairs are shown in Fig. 10. The couplings in (14) lead to the differential cross section where the Mandelstam variables are bound to satisfy s + t + u = 2(m 2 S + m 2 h ), Γ h is the total Higgs decay width and the angle θ is related to t via the relation The overall combinatorial factor 1/2 accounts for two identical final state particles. If we take the nonrelativistic limit s → 4m 2 S , which is the leading contribution to DM freeze-out, and we stay away from the Higgs resonance at m S m h /2, we find (ii) SS → ∆∆ * : Annihilations to scalar leptoquarks are tree-level processes with diagrams shown in Fig. 11. The total cross section for this process reads where N c = 3 denotes the number of colors, and 2T + 1 accounts for the LQ weak-isospin multiplicity, i.e. T = 0 for a weak singlet, T = 1/2 for a doublet and T = 1 for a triplet LQ. Away from the Higgs pole and in the non-relativistic limit we have the expression Note that this is the only annihilation process that depends on λ 5 coupling at tree-level. (iii) SS → V 1 V 2 : DM can also annihilate to SM gauge bosons, and these processes can proceed either via an s-channel Higgs mediated contribution followed by a tree (loop) mediated hV 1 V 2 vertex for massive (massless) vector bosons, or via direct loop diagrams. The associated Feynman diagrams are depicted in Fig. 12. We parameterize the SS → V 1 V 2 loop amplitude as where D V1V2 (s) are form-factors, and µ V (p, A) denotes the V -boson polarization with momentum p and index A. For gluons, A = 1, . . . , 8 and T AB = δ AB /2, while for the electroweak bosons one should replace T AB by 1. The factor δ V1V2 accounts for identical particles in the final state, being δ V V = 1 for V ∈ {γ, g, Z} and δ V1V2 = 0 otherwise. The above gauge-invariant form is valid in the m V1,2 m ∆ limit 4 , a good approximation supported by the lower bounds on m ∆ determined by direct searches for LQs at the LHC [68]. The general expressions we have obtained for D V1V2 are also reported in App. 2. By using the expressions defined above we can express the cross section for SS → V 1 V 2 in terms of D V1V2 and of the Higgs coupling to vector bosons. We start by considering the SS scattering into W W and ZZ. In this case, the Higgs-mediated diagram in Fig. 12 appears already at tree-level where V = W, Z, and δ V V is such that δ ZZ = 1 and δ W + W − = 0. The leptoquark-loop contributions D SD V V can be found in App. 2. For the remaining annihilation channels, namely SS → γγ, SS → gg and SS → γZ, the Higgs-mediated contribution appears only at one-loop level and is included in the definition of the D V1V2 form-factors. We obtained with explicit expressions for D γγ , D gg , and D γZ reported in App. 2. Figure 12. Tree-and one-loop diagrams for SS → V1V2 process, with V1,2 a SM gauge boson. The gray blob in the first diagram appears at tree-level for SS → W W and SS → ZZ, while the same couplings arises at one-loop level SS → γγ and SS → γZ. 4 In the effective theory limit, √ s, m V 1,2 m ∆ , in such a way that the above amplitude corresponds to the effective Lagrangian (iv) SS → ff : Finally, we derive the expression for the DM annihilation into fermions. Previous studies have only considered the tree-level Higgs-mediated contributions for these processes [21]. Loopinduced ones can also be relevant for f = due to a chiral enhancement (∝ m q /m ) which overcomes the suppression by the small lepton Yukawas of the tree-level diagrams. The relevant diagrams for this process are depicted in Fig. 13, which can proceed either via λ 4 S 2 |H| 2 interaction followed by the h + − vertex (left panel), or via λ 3 S 2 |∆| 2 (right panel). Figure 13. Leading contributions to the process SS → + − . In the first diagram, the gray blob also includes the LQ loop modification of the h + − coupling. See text for details.
First, we discuss the contributions from the left diagram in Fig. 13. These contributions amount to an effective modification of Higgs Yukawa coupling, where we have kept only the dominant contributions, which arise from top-quark loops, in the leading logarithm approximation. This contribution can be combined with the LQ loops depicted in the right diagram of Fig. 13, which also induce chirality-enhanced contributions, in such a way that the total cross section reads where we assume m √ s and m m q , and we separate the real and imaginary parts of the Yukawa couplings which generate non-interfering scalar and pseudoscalar amplitudes, respectively. The loopfunction G(y, x q ), with y = s/(4m 2 ∆ ) and x q ≡ m 2 q /m 2 ∆ , is where C 0 stands for the three-point Passarino-Veltman function [110], with the same conventions as used in the Package-X documentation [111]. Finally, note that in the case of SS → qq the loop contributions are expected to be less important since the chirality factor would now be m /m q . Therefore we approximate the cross section by its the tree-level contribution via the Higgs-portal (A15)

LQ corrections to Lepton Yukawa couplings
We begin with an estimate of the LQ contributions to the Higgs Yukawa coupling to leptons at one-loop. To this purpose, we integrate-out the LQs at tree-level and incorporate the leading-logarithm contribution from electroweak running [102] to the operator The presence of this operator at the electroweak scale breaks the SM linear relation between Higgs Yukawas and lepton masses, thus inducing modification of the DM annihilation cross sections. Potentially large contributions with non-chiral couplings are only from R 2 and S 1 leptoquarks which generate at tree-level a non-chiral operator, Closing the heavy-quark loop and attaching to it Higgs line(s) causes mixing into Q eH rs operator [112][113][114], where µ EW ≈ m h denotes the electroweak scale, Λ ≈ m ∆ is the matching scale where a LQ is integrated out, and we have kept only the leading-logarithm contributions in this expression. The dim-6 operator misaligns the Yukawa couplings relative to the lepton mass matrix M , such that the Yukawa coupling reads, in the mass basis, It is clear that only the dimension-6 operator contributes to the effective Yukawa coupling modification, whereas the dimension-4 Yukawa has been absorbed in the weak-to-mass basis rotation matrices. We have neglected running of the lepton masses below scale µ EW . If the effective coefficient C lequ has an imaginary part then we get also a pseudoscalar-type Yukawa coupling¯ iγ 5 h. There are two non-chiral LQ models that contribute to the C lequ (Λ) coefficients: R 2 which is a F = 0 LQ, and S 1 with |F | = 2. Their tree-level matching relations are where the matrices y L(R) are defined in Table I. Since we consider only couplings to top quarks and leptons with equal flavour, the only nontrivial flavour combination is prst = tt in both cases. Inserting the above expressions into y eff yields Eq. (A12). The non-logarithmic contributions to this matching have been recently computed in Ref. [108,115], see also [116].

SS → V1V2 and h → V1V2 form factors
We provide the form factors D V1V2 (s) that appear in the DM annihilation cross sections reported in App. A. The short-distance (SD) LQ loop contributions are represented by the last three diagrams in Fig. 12: where x ∆ ≡ s/4m 2 ∆ . In these expressions we have set the final state masses m V1,2 = 0 consistently with the gauge-invariant expression (A7). The SU (2) L Dynkin index L 2 (∆) takes values 0, 1/2, 2 in case when weak isospin of ∆ is T = 0, 1/2, 1, respectively. Another consistency check that we find valid is equality D SD γZ = 2D SD γγ in the limit s W → −1/ √ 2, c W → 1/ √ 2, g 2 → 0, and T 3 → 0, where both Z and A are massless and have identical term in the covariant derivative (1).
For SS → ZZ and SS → W W there is an additional s-channel Higgs mediated diagram which interferes with the leptoquark loop contribution D SD V1V2 , as explicitly shown in Eq. (A8). For the form factors involving a massless boson we can absorb such s-channel Higgs contributions into D V1V2 : .
The above amplitude corresponds to effective Lagrangian L hV1V2 = −(1/2)g hV1V2 hV 1µν V µν 2 . The Higgs to diboson decay widths are then The width of h → γγ has a factor of 2 relative to h → γZ due to identical particles, whereas h → gg has additional factor of 2 that stems from the gluon octet sum, A,B (δ AB /2) 2 = (N 2 c − 1)/4. The couplings g hV1V2 in the SM have been computed in [117][118][119] and consist of loop contributions dominantly from gauge bosons and top quark. The leptoquark contributions are analogous to the diagrams of SS → V 1 V 2 process shown in the rightmost three diagrams in Fig. 12, where we have to replace SS∆∆ * with the h∆∆ * vertex. This similarity between the two processes allows us to identify g hV1V2 = − 2λ5v λ3 D SD V1V2 . The expressions for g hV1V2 are where x i = s/(4m 2 i ) and λ i = (4m 2 i )/m 2 Z , for i ∈ {W, t, ∆}, and the loop-functions are reported in App. 3. In these expressions, we have neglected the sub-leading contributions from light fermion loops. In these expressions we can further simplify sums over the weak isospin states of ∆: and, similarly,