B anomalies and dark matter: a complex connection

We study an extension of the Standard Model that addresses the hints of lepton flavour universality violation observed in $B\to K^{(*)} l^+l^-$ decays at LHCb, while providing a viable candidate for dark matter. The model incorporates two new scalar fields and a Majorana fermion that induce one-loop contributions to $B$ meson decays. We show that agreement with observational data requires the new couplings to be complex and that the Majorana fermion can reproduce the observed dark matter relic density. This combination of cosmological and flavour constraints sets an upper limit on the dark matter and mediator masses. We have studied LHC dijet and dilepton searches, finding that they rule out large regions of parameter space by setting lower bounds on the dark matter and mediator masses. In particular, dilepton bounds are much more constraining in a future high-luminosity phase. Finally, we have computed the scattering cross section of dark matter off nuclei and compared it to the sensitivity of current and future direct detection experiments, showing that parts of the parameter space could be accessible in the future to multi-ton experiments. Future collider and direct DM searches complement each other to probe large areas of the parameter space of this model.


Introduction
LHCb has reported anomalies in the measured decay rates of the B meson, which have been interpreted as hints of lepton flavour universality violation [1,2]. The SM predicts equal rates for the processes B → K ( * ) µ + µ − and B → K ( * ) e + e − , and it is customary to study the ratios of these branching ratios, defined as R(K) and R(K * ), since the dependencies on hadronic matrix elements (and associated uncertainties) cancel out [3]. The measurements of these hadronically clean observables deviate consistently (although perhaps with not enough statistical significance) from the SM prediction R(K ( * ) ) = 1 [4]. These hints are complemented by measurements of other observables that are more sensitive to hadronic physics. In particular, the differential branching fractions [1,2,5] and angular observables [6][7][8][9][10][11][12][13] associated to the processes B → φµ + µ − and B → K ( * ) µ + µ − also deviate from the SM predictions. Interestingly, all the apparent anomalies involve the transition b → sµ + µ − .
In order to account for these experimental results, one can modify the SM effective Hamiltonian, which involves penguin and box diagrams, by including one-loop contributions from new exotic particles. A full classification of the various particle combinations, considering different gauge representations, was presented in Refs. [14,15]. Among the different models, some featured neutral scalar or fermions that, if stable, could play the role of dark matter (DM) 1 . The first possibility was investigated in Ref. [41], where it was found that the large new couplings required to reproduce the correct DM relic abundance induce sizeable 1-loop contributions to DM-nucleon scattering, leading to very strong limits from direct detection experiments. In addition, as reported by [42], the Higgs portal coupling typically dominates over other new physics effects. The second possibility was addressed in Ref. [43], where the fermionic dark matter field was accompanied by one additional scalar and one additional coloured fermion.
In this work, we consider a modification of the model of Ref. [43]. Namely, we will also assume a fermionic dark matter particle, but with two extra scalar fields, one of which has a colour charge. On top of this, we include the latest SM theoretical prediction for the mass difference in B s −mixing [44], which differs from the experimental observation by 1.8 σ.
In order to reduce this tension and provide an explanation for the B anomalies, complex couplings are needed, leading to new CP-violation sources, a scenario that has not been studied in the context of one-loop models so far. We explore the parameter space of this model, taking into account all the flavour observables, DM constraints, and LHC collider signatures. This paper is organised as follows. In Section 2, we introduce the details of the particle physics model, address the constraints from the observed DM relic abundance and B s −mixing and discuss the implications on the model's parameter space. In Section 3, we investigate the possibility of observing this scenario at the LHC, for which we take into account dijet and dimuon searches. We also include a projection of the potential reach of the High Luminosity phase of the LHC. Finally, in Section 4, we compute the DM-nucleus scattering cross section and study current constraints and the future reach of direct DM detection experiments. The conclusions are presented in Section 5.

The model
In this article, we consider a model in which the DM particle is a Majorana fermion, χ, with two extra scalar fields, φ q and φ l , which couple to left-handed quarks and leptons, respectively 2 . The interactions between the new particles and the SM are described by the Lagrangian, where Q i and L i denote the SM left-handed quark and lepton doublets of each generation, and λ Q i and λ L i are the corresponding new couplings. The quantum numbers for the new fields are summarised in Table 1. We impose a Z 2 parity under which the SM fields are invariant, and which guarantees the stability of the DM candidate, as long as m φ q,l > m χ . Upon rotation from the electroweak to the quark mass eigenbasis, the couplings λ Q i are rotated in flavour space. Assuming that the electroweak and mass eigenbasis are aligned for the leptons and down-type quarks, the couplings to the up-type quarks are generated by the CKM rotation as follows: From now on, we will denote the couplings in the mass eigenbasis with the corresponding quark or lepton label. These couplings are, in general, complex. This model induces new physics contribution to flavour observables at the one loop level. In particular, a new box diagram appears for the b → sµ + µ − transition, as shown in Figure 1. These effects can be described using an effective field theory approach, thus parameterising the new contributions as corrections to the corresponding Wilson coefficients of the effective Hamiltonian, The Wilson coefficients C 9 , C 9 , C 10 , C 10 contain both the SM and new physics (NP) contributions, with the primed coefficients defined in an equivalent way.
Global fits [20,[45][46][47][48][49][50][51][52][53] have been used to determine the new physics contribution to the Wilson coefficients in order to reproduce the observed experimental results. These fits favour C NP 9 = −C NP 10 , and suggest that no new physics is required for operators involving electrons or tau leptons. Motivated by these results, we assume negligible couplings to the first quark generation (i.e., λ Q 1 = 0) and to the first and third lepton generations (i.e., λ e = λ τ = 0). This provides an explanation for the R K ( * ) anomalies, while relaxing the bounds from other searches. Therefore, in total, we are left with six free parameters in this model, namely the masses of the three new particles (m χ , m φ l , m φq ), and the couplings to b−type quarks, s−type quarks, and leptons (λ b , λ s , λ µ ).
It should be noted that the couplings λ 1 |φ l | 2 |H| 2 and λ 2 |φ q | 2 |H| 2 are allowed by gauge symmetry in the Lagrangian of Eq. (2.1). However, they only lead to an overall shift to the masses of φ l and φ q after electroweak symmetry breaking since the couplings to the Higgs play no phenomenological role in the relevant range of φ l,q masses. Likewise, the terms λ 3 |φ l H| 2 and λ 4 |φ q H| 2 are also allowed by gauge symmetry. They typically induce a small split in the masses of the neutral and charged components of the doublets φ l and φ q in the range of φ l,q masses that survive the collider constraints. Finally, a term of the form (φ l H) 2 can lead to large contributions to neutrino masses at one loop, which forces the corresponding coupling to be extremely small [43]. We will neglect these couplings in the following.
As mentioned in the Introduction, similar models have been discussed in the literature, featuring either scalar DM [41,[54][55][56] or fermionic DM [43]. Our model differs from that of Ref. [43] in that we have two extra scalar fields which couple to the lepton or quark sectors.

Dark matter relic abundance
In order for χ to be a viable DM candidate, it must reproduce the observed relic abundance, which can be inferred from Planck satellite data to be Ωh 2 = 0.1199 ± 0.0022 [57]. The pair-annihilation proceeds through the two t−channel diagrams with φ q and φ l , shown in Figure 2.
The stringent flavour constraints force the couplings to quarks to be much smaller than the couplings to leptons (muons and neutrinos), and the combination of flavour and collider bounds impose m φq > m φ l , with coloured scalars generally above 1 TeV. Therefore DM annihilation into a µ − µ + or ν µνµ pair is the dominant channel. The thermally-averaged annihilation cross section, σv , can be expressed as a plane wave expansion in terms of the dimensionless parameter x = m χ /T . For the case of a Majorana fermion, the zerovelocity term is helicity suppressed, and the leading contribution comes from the linear term in 1/x [58], where we have neglected the muon and the neutrino masses. In order to reproduce the correct relic abundance, we can now impose σv = 2.2 × 10 −26 cm 3 s −1 (where x ∼ 20 at freeze-out).
We will use this relation to fix m φ l as a function of the other parameters, thus effectively reducing by one the number of free parameters. Furthermore, due to the suppression of the velocity-independent term for σv , indirect detection bounds are not expected to constrain our model.

B s −mixing and other flavour constraints
This model introduces new couplings to the s and b quarks (and to the rest of the quarks by rotation of the CKM matrix). We must therefore incorporate constraints from B meson physics.
The most relevant bounds are those that involve b → sµ + µ − transitions. The new physics contribution to the Wilson coefficient comes from box and photon-penguin diagrams [14,15], where we have defined the dimensionless variables x q = m 2 φq /m 2 χ and x l = m 2 φ l /m 2 χ , and the loop functions are: , The term G(x q , x l ) vanishes if χ is a Dirac particle.
In order to constrain the Wilson coefficients we use the first global fit that takes into account the possibility that C 9 and C 10 are complex [59]. This is a natural scenario that arises when new CP-violation sources are introduced, and has not been studied in detail in the literature so far.
Likewise, the new physics contribution to B s −mixing can be parameterised in terms of an effective Hamiltonian, where α and β are colour indices. The new physics contribution to the Wilson coefficient is given by where the loop functions F and G were already defined in Eq. (2.11).
In order to quantify the allowed magnitude of the Wilson coefficient C NP BB , we follow the steps of [44] and introduce a complex parameter ∆ in the following way: where M SM 12 and M NP 12 describe the SM and new physics contributions to B s −mixing, and their values are given by the corresponding box diagrams. The complex phase, φ ∆ , quantifies the CP-violating effects introduced by the imaginary parts of the new couplings. We find: where ∆M s is the mass difference of the mass eigenstates of the B s meson.
The parameter |∆| can be constrained using the most precise experimental measurement of ∆M s [60] and the last update on its theoretical prediction [44], which show a 1.8σ difference, (2.16) The dominant uncertainties in the calculation of ∆M SM s come from lattice predictions for the non-perturbative bag parameter, B, and decay constant, f Bs , and to a lesser extent from the uncertainty in the values of CKM elements. Both of these errors have been considerably reduced since the last theory update for the mass difference [61]. The last average given by the lattice community [62] gives significantly more precise values for B and f Bs .
From these values, one can infer |∆| = 0.887 ± 0.055, and using the data provided in Ref. [44] we obtain C SM BB = 4.897 × 10 −5 TeV −2 . Using Eq. (2.15) we find that the Wilson coefficient has to satisfy CP-violating effects are further constrained by the CP asymmetry of the golden where β s = 0.01852±0.00032 [63], and penguin contributions are neglected. Using Eq. (2.15), this can be interpreted as an additional constraint on the real and imaginary parts of C NP BB (and in turn, on the real and imaginary parts of the couplings λ s λ * b ).
In Figure 3, the effect of all of these constraints on the real and imaginary parts of the couplings λ s λ * b for several benchmark points is shown. Regions that are allowed by b → sµ + µ − observables and B s −mixing (given by Eqs. (2.17) and (2.18)) are shaded in green and blue, respectively. For illustrative purposes, the figure shows the constraints for multiple values of DM and mediator masses, while keeping λ µ = √ 4π fixed. We remind the reader that the mass of m φ l is fixed so as to reproduce the correct relic density using Eq. (2.9).
As we can observe, in order to simultaneously satisfy both types of constraints, complex couplings are needed (Im(λ s λ * b ) = 0). Also, as the mass of the dark matter particle and the mediators increase, both areas are more difficult to reconcile. In practise, this leads to an upper bound on the masses of the exotic new particles. The precise limit depends on the choice of couplings, which we will discuss in Section 3.
Finally, the new physics couplings to the up-type quarks are generated via CKM rotation, These couplings generate a new physics contribution to D 0 −mixing, and the Wilson coefficient C NP DD is obtained replacing λ s and λ * b in Eq. (2.13) by λ u and λ * c , respectively.
In contrast to B s −mixing, there is no precise theory determination for the mass difference in the D 0 system. Therefore, in order to constrain the new physics contribution to C DD we whereas the new physics contribution to D 0 −mixing is described by where O is a combination of operators containing all possible SM and new physics contributions to D 0 −mixing. Using the last results from [65] we get the following bound on the Wilson coefficient: Although this model induces new physics contributions to other flavour observables (such as b → sγ, b → sνν and effective Zµ + µ − and Zq i q j couplings), their size is very small and does not produce significant deviations from current experimental searches.

Benchmark scenarios
All the new physics contributions to the observables described above depend on five independent parameters: the three masses of the new particles, m χ , m φq and m φ l , the product of the couplings λ s λ * b and the absolute value of the coupling |λ µ |.
The three masses only enter the Wilson coefficients through the factor m −2 χ and the dimensionless loop functions. In addition, all the Wilson coefficients are proportional to λ s λ * b or |λ µ | 2 or both. In order to constrain our model, we consider two scenarios by fixing the value of |λ µ |. Then we scan over the mass parameters m χ and m φq , with m φ l fixed by the requirement of reproducing the correct relic abundance, and check all the flavour observables described in Section 2.2. In this way, for any combination of masses and a fixed value of |λ µ | we get a set of allowed values for λ s λ * b . We consider two hierarchies between |λ s | and |λ b | that lead to different constraints from D 0 −mixing, and, ensuring that Im(λ s λ * b ) = 0, we define the following benchmark scenarios: where |λ µ | = √ 4π is the perturbative limit. After establishing a hierarchy between |λ s | and |λ b |, we calculate their maximum and minimum allowed values from the corresponding maximum and minimum allowed values of λ s λ * b . Scenarios with |λ s | > |λ b | are excluded by D 0 −mixing constraints. Likewise, as we will see in Section 3, smaller values of λ µ are constrained by LHC bounds.

LHC constraints and prospects for high-luminosity
In this section, we study the experimental signatures that this model would produce at the LHC. DM search strategies in both ATLAS and CMS involve analysing final states containing jets and leptons produced in association with a DM particle, identified from missing transverse energy. In this model, direct production of the coloured and leptonic scalar doublets φ q and φ l , respectively, typically leads to such final states.
Let us first consider production processes that involve the coloured scalar, φ q . In this case, our model could lead to visible signals in final states with both monojet / dijet + E T / signatures. When the new physics coupling λ q is smaller than the strong interaction coupling, α QCD , pure QCD processes constitute the main contribution to the cross section [66]. In this model, this implies that QCD diagrams dominate over those with new physics couplings. As a consequence, monojet searches for this model are less effective than dijet searches and we will concentrate on the latter. The dijet + E T / processes are shown in Figure 4, where diagrams (a) correspond to the QCD contributions, and diagrams (b) and (c) involve new physics couplings. The main production channel is the pair production of the coloured scalar particles, that subsequently decays into a DM particle and a quark, In addition, the scalar doublet φ q has the same quantum numbers as squarks in supersymmetric (SUSY) models. Therefore, the kinematics in its production and decay in diagrams (a) of Figure 4 mimic those of squarks in SUSY models with decoupled gluinos. As a consequence, limits from ATLAS and CMS squark searches can be used to constrain the model.
One can also consider the pair production of the leptonic scalar, φ l . In this case, the production process is mediated by W or Z bosons and involves the electroweak coupling, as shown in Figure 5. The decays of φ l lead to clean final states with one or two leptons and missing energy. Although flavour constraints require λ µ λ q , the cross section of this process is smaller than the production of the coloured mediator for similar mediator masses. However, since m φ l is fixed for every value of m χ to reproduce the correct relic abundance, there are regions of the parameter space where both searches are complementary. We will here consider the process where the dimuon channel leads to the strongest constraints. As in the previous case, we can exploit the analogy between φ l and sleptons to use the limits from slepton searches to constrain this model.

Simulation details
We have implemented this model in Feynrules 2.3 [67]. The calculation of the matrix elements and the event generation is done using MadGraph5 aMC@NLO 2.6.3 [68]. Production and decay of the new particles are considered independently using the narrow width approximation, as implemented in MadSpin [69], which further accounts for spin correlations in decay chains 4 . We then use Pythia 8.235 [70] to shower the parton-level events and we pass the output to CheckMATE 2.0.26 [71], which compares the expected signal with supersymmetric searches at the LHC and derives an exclusion limit. As we have explained above, we can apply squark and slepton searches to constrain the coloured and leptonic mediator, respectively.
In order to describe initial and final state radiation and reproduce the correct jet struc-ture precisely, we consider leading order (LO) production with parton shower matching and multijet merging when needed. The LO multijet merging techniques describe how parton shower emissions can be combined with full matrix element calculations to achieve a better accuracy in the description of the radiation spectrum. Using this technique, every jet is classified according to its p T and then compared to a hardness scale Q cut . In this way, emissions above the hardness scale Q cut are described at LO accuracy using the corresponding matrix element calculation for an extra hard, wide-angle QCD emission in the final state, while emissions below this scale are defined as soft or collinear jets and the all-orders resummation description from the parton shower is preserved. Note that even though O(α s ) corrections are included using this procedure, the calculation remains formally LO + LL accurate after parton shower due to missing virtual corrections.
After hadronization, the showered events and the production cross sections are passed to CheckMATE. Each model point is tested against all the implemented experimental analyses to determine the optimal signal region. For this signal region, CheckMATE compares the simulated signal with the actual experimental observation and determines whether the model point is excluded at the 90% confidence level.

Results
Constraints from LHC searches for the four benchmark points defined in Section 2.3 are presented in Figure 6 on the (m χ , m φq ) plane, for all the points that satisfy the flavour constraints of Section 2.2 and that reproduce the correct DM relic abundance. This figure shows the complementarity between the experimental limits obtained from the pp → jj + E T / and pp → µµ + E T / searches. The experimental results used in our analysis are summarised in Table 2. The colour code represents the average value of the coupling |λ b | in the region allowed by flavour constraints, defined as |λ b | mean = (|λ b | max + |λ b | min )/2, where |λ b | max and |λ b | min are the maximum and minimum allowed values respectively. The variation of our results when choosing either the minimum or maximum value for |λ b | has been checked and is insignificant.
Regarding the pp → jj + E T / search, the limits in every scenario show that for the lightest DM mass, coloured mediators with masses below ∼1 TeV are excluded. Even though heavier DM produces larger amounts of missing energy in final states, the cross section decreases rapidly with the m χ , leading to similar exclusion limits. It is interesting to note that exclusion limits are slightly stronger for the scenarios with |λ b,t | > |λ s,c |, where mediators with masses below ∼1.1 TeV are excluded. The reason for this is that final states with either top or bottom quarks are more sensitive to some experimental searches. The most stringent experimental search involves final states with at least two (bb production) or four (tt production) jets or exactly two leptons and missing energy [72]. In particular, the most sensitive signal region is optimised to detect events featuring a DM particle produced in association with a tt pair, which decays fully hadronically.
Regarding the pp → µµ + E T / search, the limits show that models with dark matter masses below approximately 30 GeV are ruled out for |λ µ | = 2, with the exclusion limit going down to ∼ 13 GeV for |λ µ | = √ 4π. This corresponds to mediator masses below 360 GeV for |λ µ | = 2 and 410 GeV for |λ µ | = √ 4π. The pp → µµ + E T / cross section mainly depends on m φ l , so the limits on m χ can be understood through its relation with m φ l given by the DM relic condition (2.9) for a particular value of λ µ . The most stringent search involves final states with 2l + 0j, 2l and at least 2 jets, or 3l and missing energy [73]. In particular, the most sensitive signal region is characterised by 2l + 0j and a dilepton invariant mass 13 TeV 1710.11412 [72] tt / bb + E T / 36.1 ATLAS-CONF-2017-039 [73] 2l / 3l + E T / 36.1 14 TeV atlas phys pub 013 011 [74] tt + E T / 3000 atlas phys 2014 010 hl 3l [75] 2l / 3l + E T / 3000 Table 2: List of experimental searches sensitive to our model, where l denotes electron and muon. The third column describes the final state targeted by the analysis and the last column displays the total integrated luminosity.
m ll > 300 GeV, and it is optimised to target slepton pair production.
The most remarkable result is that LHC limits completely exclude the scenario with |λ µ | = 2 and λ b = λ * s , as well as a sizeable region of the scenario with λ b = 4λ * s for the same |λ µ |. These constraints become weaker for larger values of |λ µ | and, for the scenarios with |λ µ | = √ 4π, most of the parameter space is allowed. It is crucial to note that the limits coming from final states with jets and leptons are complementary to each other. While the former exclude regions of the parameter space with large m χ and small m φq , the latter rule out models with very heavy mediator masses m φq and light dark matter. Importantly, these limits are also complementary to the ones coming from direct detection, where dark matter masses below 12 GeV lie below the neutrino floor. Therefore, it is fundamental to consider both approaches to explore the model.
It is worth mentioning that the small couplings required by flavour constraints lead to decay widths slightly below the QCD scale for m φq 370 GeV. Strictly speaking, this means that the computation of the decay width cannot be handled perturbatively and that the new particle φ q may hadronize into bound states with SM quarks, analogous to R-hadrons [76], before decaying. However, the typical width involved is Γ φq ∼ O(10 −2 ) − O(10 −3 ) GeV, which means lifetimes of the order τ ∼ 10 −22 s, so any potential bound state would decay promptly in the detector. This region of the parameter space is excluded by ATLAS and CMS R-hadron searches [77,78].
We have also studied the limits that could be obtained with 3000 fb −1 of 14 TeV data once the LHC High Luminosity phase [79] is completed. As we can observe in the plots, the main gain would come from the leptonic channels, which would allow to test a considerable amount of the model's parameter space. In particular, scenarios with |λ µ | < 2 would be completely excluded. The experimental searches giving the strongest exclusion limits target the same final states and are shown in the low panel of Table 2.

Direct DM detection prospects
Finally, in this section we discuss whether our model is expected to produce an observable response in direct detection experiments. We have calculated this response, by matching the model parameters to effective DM-nucleon interaction terms, 1) where N is the corresponding nucleon, and O i is the set of non-relativistic operators [80,81].
The values for the coefficients c N i can be derived as the non-relativistic limit of the original interaction Lagrangian, and the differential rate can be computed using the corresponding nuclear form factors from Refs. [81,82], and for a given choice of the DM halo properties. We have adopted the so-called standard halo model [83] with local DM density ρ χ = 0.4 GeV/cm 3 , a central velocity of v 0 = 220 km s −1 , and a escape speed of v esc = 544 km s −1 to calculate the number of expected recoils in a specific experiment.
The leading tree-level DM-quark interactions are given by scalar (χχψψ) and vector (χγ µ χψγ µ ψ) type interactions. The latter is the leading contribution to O 1 for Dirac DM [84], but it vanishes in the case of Majorana DM. For scalar type interactions Majorana DM does not in general vanish, but with our models chiral structure, it does. With subdominant couplings to the first generation of quarks, and given that m φq > m φ l , one-loop contributions to the DM-nucleon scattering cross sections will generally be larger than the tree level process. The loop contributions for a generic fermionic DM that involve the exchange of a photon can be classified as electric and magnetic dipoles (χiσ µν γ 5 χF µν and χσ µν χF µν , respectively), anapole (χγ µ γ 5 χ∂ ν F µν ), and charge radius (χγ µ χ∂ ν F µν ). However, in the particular case of Majorana DM considered in this work, the magnetic dipole and charge radius effective couplings are forbidden by charge conjugation symmetry. Thus, the dominant one-loop interaction is the anapole moment [85]. When taking the non-relativistic limit, the anapole moment gives contributions to the O 8 and O 9 operators [86,87], which are velocity and momentum dependent. In terms of the fundamental parameters of the model, the corresponding couplings read where e is the electron charge, Q N is the nucleon charge, and g N are the nucleon g-factors (g p = 5.59 and g n = 3.83). The effective coupling to the anapole interaction term, A, reads [88] A = − e |λ µ | 2 96π 2 m 2 with µ ≡ m 2 φ l /m 2 χ and ≡ m 2 l /m 2 χ . The nuclear responses to the O 8 and O 9 operators are markedly weaker than that of O 1 , which implies that, in general, the scattering cross section is very small and beyond current experimental limits. Furthermore, because our DM particle interacts with the quark sector, it is not a priori clear that the spin-independent O 1 and spin-dependent O 4 arising from the so-called twist-2 operator [89][90][91] and the axial vector operator respectively are still negligible.
Given the range of DM masses that we consider in this study, the main constraint is due to Xenon1T results [92], which we simulate using the prescription outlined in appendix A of Ref. [93], achieving good agreement. As we can see in Figure 7, the theoretical predictions for this model are beyond the reach of current experimental searches. We also show the reach of future direct detection experiments. The LZ detector, will employ 5.6 tons of liquid xenon with 1000 days exposure as outlined in [94,95]. The DarkSide-20k experiment [96], is an argon detector which will employ 20 tons of fiducial mass for a duration of 10 years. We have assumed that the DarkSide collaboration will be able to achieve a threshold energy of 5 keV, a reasonable assumption considering the results from DarkSide-50 [97]. For reference we have also calculated the neutrino floor for anapole interactions in the (A, m χ ) plane. We have used the prescription described in Ref. [98] and the expected neutrino fluxes from Refs. [99][100][101][102][103][104]. It is clear that our model favourably lays in a region of parameter space that would be probed by a generation of experiments with multi-ton targets, that can probe near or even slightly beyond the neutrino floor. Spectral analysis with the neutrino background compounded with annual modulation data, could provide complete discrimination between model and the anapole moment which is both velocity and momentum dependent.
For completeness, we have also calculated the effect on the total scattering cross section from aforementioned twist-2 operator and spin-dependent interaction. The former contribute to the spin-independent scattering cross section (operator O 1 ) and can be sizeable if the new coupling to quarks is large or the colour mediator is very light. We have explicitly checked that once LHC constraints are included in the parameter space of the model, these terms are always subdominant to the anapole term discussed above. We represent in Figure 8 the theoretical predictions for c 1 as a function of the DM mass from this contribution. For the spin-dependent interaction, we found that the predicted rate for our sampled parameter space is always sub-dominant.
Had we chosen to work with a Dirac fermion, the dipole and charge radius contributions should have been added. As it has been pointed out in Ref. [84], the fairly large coupling to muons that is required to explain the flavour anomalies leads to effective DM couplings that are orders of magnitude higher than those coming from the tree level contribution, the most important being the charge-radius interaction. This we have checked, and in fact above m χ ∼ 10 GeV, all our parameter points are excluded by Xenon1T. Below m χ ∼ 10, the model For comparison, we show the current exclusion line by Xenon1T [92] and the predicted reach of LZ [94,95] and DarkSide-20k [96]. The shaded area represents the neutrino floor. The plot on the right-hand side incorporates LHC constraints, explained in more detail in Section 3.  Figure 7, but for the spin-independent coupling, c 1 , that originates from the twist-2 coupling. is excluded by both LHC constraints and indirect detection bounds. Unlike in the Majorana case, the s-wave contribution to the thermal cross section σv is no longer helicity suppressed and hence excluded [105].
Our results suggest that future multi-ton direct detection experiments, such as DarkSide [96], would be able to probe this model in the mass range m χ ∼ 10 − 60 GeV. It is very interesting to point out that many of the points in this DM mass range feature very heavy φ q and therefore would be beyond the reach of collider searches. In a sense, future direct DM detection and the LHC complement each other to probe a large part of the model's parameter space.

Conclusions
In this article, we have studied a particle physics model that addresses the hints of lepton flavour universality violation observed by LHCb in b → sµ + µ − transitions, and that provides a solution to the dark matter problem. The scenario that we have analysed incorporates two new scalar fields and a Majorana fermion that provide one-loop contributions to B meson decays.
The Majorana fermion is stable and can reproduce the observed DM relic abundance. We have studied the effect of new physics in flavour observables, for which B s −mixing and b → sµ + µ − processes provide the most important constraints. In order to find an explanation for the B anomalies and to reduce the 1.8 σ tension between the predicted and measured mass difference in B s −mixing, complex couplings are needed. We have used results from the first global fit that takes into account this possibility. The combination of flavour bounds and constraints on the DM relic abundance leads to upper limits on the masses of the exotic states, and in general points towards a rather light DM candidate (with a mass m χ 200 GeV).
We have studied the signatures that this model would produce at the LHC. The dominant processes are the pair production of the coloured and leptonic scalars. For the former, the strongest exclusion limits are given by dijet + E T / searches. For the latter, the final states are very clean, containing 1 or 2 leptons and missing energy. Both searches are complementary and exclude different regions of the parameter space, setting lower bounds on DM and mediator masses. The high-luminosity phase improves bounds coming from both searches, with dilepton being the most pronounced. The collider constraints are weakened when the λ µ parameter is pushed towards the perturbative limit.
Finally, we have investigated how DM direct detection experiments constrain this model. Given the range of DM masses that we consider in this study, the main constraint is due to Xenon1T results. The small new couplings required by flavour constraints means that one-loop contributions to the DM-nucleon scattering cross section are generally larger than the tree level process. In particular, the dominant loop induced interaction is the anapole moment. We have shown that this model is not excluded by current data and could be probed by the next generation of experiments with multi-ton targets in the mass range m χ ∼ 10 − 60 GeV.