Systematic approach to $B$-physics anomalies and $t$-channel dark matter

We study renormalisable models with minimal field content that can provide a viable Dark Matter candidate through the standard freeze-out paradigm and, simultaneously, accommodate the observed anomalies in semileptonic $B$-meson decays at one loop. Following the hypothesis of minimality, this outcome can be achieved by extending the particle spectrum of the Standard Model either with one vector-like fermion and two scalars or two vector-like fermions and one scalar. The Dark Matter annihilations are mediated by $t$-channel exchange of other new particles contributing to the $B$-anomalies, thus resulting in a correlation between flavour observables and Dark Matter abundance. Again based on minimality, we assume the new states to couple only with left-handed muons and second and third generation quarks. Besides an ad hoc symmetry needed to stabilise the Dark Matter, the interactions of the new states are dictated only by gauge invariance. We present here for the first time a systematic classification of the possible models of this kind, according to the quantum numbers of the new fields under the Standard Model gauge group. Within this general setup we identify a group of representative models that we systematically study, applying the most updated constraints from flavour observables, dedicated Dark Matter experiments, and LHC searches of leptons and/or jets and missing energy, and of disappearing charged tracks.


Introduction
The first decade of operation of the Large Hadron Collider (LHC) has resulted in the tremendous success represented by the discovery of the Higgs boson and provided us with a host of precise measurements and searches for new phenomena, finding no conclusive evidence of departures from the predictions of the Standard Model (SM). Nevertheless, the SM leaves unanswered a number of fundamental questions that provide strong motivation for its extension. The most compelling problem is that the SM lacks a candidate of dark matter (DM), whose existence has been established by an impressive number of cosmological and astrophysical observations, spanning many orders of magnitude in redshift: from the Cosmic Microwave Background (CMB) to galactic rotation curves, see Ref. [1] for a review. However, we do not have at the moment knowledge about the nature of dark matter nor about its mass and interactions with the SM sector, being all evidence based on its gravitational effects. All direct and indirect searches for particle dark matter have so far given negative results. Nevertheless, it is plausible that DM interacts to some extent to the SM fields, as a substantial DM abundance must be produced in the early universe. This is the case of the thermal freeze-out mechanism, which assumes that DM is a thermal relic, most commonly a weakly-interacting massive particle (WIMP). Sizeable interactions to SM particles are then required to keep DM in thermal equilibrium with the SM bath in the early universe and to ensure an efficient DM annihilation mechanism in order to avoid the WIMP relic density to be larger than the DM abundance that Class F − Fermion mediator Class S − Scalar mediator is observed today. In this work, we are going to assume that the observed DM density is accounted for by a thermal WIMP that interacts with SM quarks and leptons and other extra fields in a way that can address the so-called B-physics (or flavour) anomalies. In fact, although direct searches performed by the LHC collaborations for the production of new particles have found no evidence of their existence, several experimental collaborations, with LHCb being the prominent one, have brought to light a persistent and coherent pattern of deviations from the SM predictions in semileptonic decays of B mesons of the kind b → s + − . This could very well be the first experimental hint for beyond the SM (BSM) physics at energies not much larger than the electroweak scale. In particular, LHCb and B-factory experiments observe a deviation from Lepton Flavour Universality (LFU) predicted by the SM in the theoretically clean observables R K ( * ) ≡ BR B → K ( * ) µ + µ − /BR B → K ( * ) e + e − [2][3][4]. Moreover, a number of measurements are in tension with the SM predictions for the branching ratios and angular distributions of several b → sµ + µ − modes [5][6][7][8][9][10][11][12][13][14]. All these anomalies could be explained by a deficit of b → sµ + µ − events compared to SM expectations due to the interference between SM and BSM amplitudes. The simplest way to achieve such an effect is to add non-standard contributions, δC 9,10 µ , to the following operators e 2 16π 2 V tb V * ts C 9 µ (sγ µ P L b)(µγ µ µ) + C 10 µ (sγ µ P L b)(µγ µ γ 5 µ) + h.c. .
While not providing the absolute the best fit to the anomalies, an interesting scenario, still in excellent agreement with the data, is represented by δC 9 µ = −δC 10 µ ≈ −0.5, corresponding to the case of only lefthanded (LH) currents entering Eq. (1). According to global fits to B-physics data, such a scenario is preferred to the SM prediction at the ∼ 5σ level [15][16][17][18][19][20][21][22][23][24][25][26]. This does not reflect of course an established breakdown of the SM: a combination of overlooked systematics, statistical fluctuations, and underestimated hadronic uncertainties could conspire to account for such a large deviation from the SM in the global fit. Nevertheless it is tempting to explore new physics (NP) scenarios that could explain the anomalies and assess their capability of addressing other shortcomings of the SM, in particular the DM problem.
In this paper we systematically build a set of simplified models that can explain the B-physics anomalies and simultaneously provide a good DM candidate, and we study their phenomenology with a particular focus on the LHC limits on production of new heavy particles and the bounds from direct-and indirectdetection DM searches. Our aim is to highlight the minimal building blocks that a more complete theory may need to include. For the sake of minimality we are going to employ the following procedure. Figure 2: Illustrative DM (co-)annihilation diagrams for the case of fermion DM belonging to the field denoted as Ψ in Figure 1. Analogous diagrams arise in the other cases. Gauge diagrams such as the third one are only present if DM belongs to an SU (2) L multiplet (Ψ ± are charged states in the same multiplet).
• We focus on minimal solutions of the B-physics anomalies of the δC 9 µ = −δC 10 µ kind, hence we only introduce new fields -in the lowest possible number -that couple to left-handed quarks and leptons (the SM SU (2) L doublets).
• We require that at least one of the BSM fields contains a state which can be a good DM candidate, i.e. neutral and colour singlet.
• We assume that DM stability is ensured by an unbroken symmetry (a Z 2 parity or another global symmetry), which forbids interactions between a NP particle and two SM particles, as well as mixing between NP and SM fields. As a consequence the BSM contributions to C 9 µ and C 10 µ can only arise through one-loop diagrams like those shown in Figure 1 with only BSM fields running in the loop, as in the framework studied in Refs. [27][28][29]. Notice that for simplicity we only consider spin 0 and spin 1/2 fields and that only three new fields need to be added to the SM.
• SU (3) c ⊗ SU (2) L ⊗ U (1) Y gauge invariance and the requirement of a consistent DM candidate tightly constrain the possible quantum numbers of the BSM fields. Furthermore, imposing the predicted relic density to be at (or below) the observed value Ω DM h 2 0.12 results in non-trivial conditions on the spectrum and couplings of the new particles, such that DM efficiently annihilates into SM particles. An unavoidable annihilation mode is given by the t-channel exchange of the other fields entering the loop of Figure 1, possibly alongside coannihilations and processes involving gauge interactions (if DM belongs to a non-trivial representation of the electroweak gauge group), see Figure 2.
To the best of our knowledge, this is the first systematic study of the connection between flavour anomalies and dark matter. However, several previous works in the literature proposed specific models that fulfill the conditions outlined above, which then will be included in our classification, see Refs. [30][31][32][33][34][35]. Other works that addressed simultaneously DM and the B-physics anomalies (among other observables) include Refs. .
The outline of the paper is the following. In Section 2 we define our setup and show the set of minimal models that follows from the approach described above. In Section 3 we describe our strategy and how we impose constraints from B-physics, LHC searches, and DM phenomenology. In Section 4 we choose a number of representative models, for which the results of our analysis are presented in detail. Finally we summarise and conclude in Section 5.

Setup
As discussed above, we consider models that can give rise to the contributions to b → sµµ processes shown in Figure 1. 1 We classify our models in two classes according to the spin of the field that couple to both quarks and leptons -which we call "flavour mediator" independently of its quantum numbers -as follows.
Class F. These models feature a vector-like fermion Ψ as flavour mediator and two extra scalars Φ Q and Φ L coupling to the SM left-handed fermions with interactions described by the following Lagrangian: Class S. In these models, we introduce a scalar flavour mediator Φ and two fermions Ψ Q and Ψ L in vector-like representations of the SM gauge group: In the spirit of our simplified-model approach, we are considering non-zero couplings of the new fields only to second and third generation left-handed quarks (Γ Q 2 , Γ Q 3 ) and muons (Γ L 2 ). For more definiteness we will use, throughout the paper, the following notation: These couplings are defined in the basis where the down-quark and charged-lepton mass matrices are flavour diagonal. Furthermore we assume a global symmetry, whose effect is to forbid mixing between extra fields and SM fields and ensure that the lightest new state is stable. This can be achieved by introducing an unbroken Z 2 parity under which the SM fields are even and the new fields are odd, or an equivalent continuous symmetry. Finally, unless Table 2: Models with a fermion flavour mediator (Class F). The fields are denoted by their transformation properties under, respectively, (SU (3) c , SU (2) L , U (1) Y ). We highlight in cyan the models that we study in detail in Section 4. otherwise stated, we will usually assume the interactions in the scalar potential -such as the unavoidable quartic couplings between our new scalars Φ X and the Higgs field H of the form Φ † X Φ X H † H -be small enough to have only subdominant effects on the phenomenology of our models.
The possible gauge quantum numbers of the extra fields follow from the requirement of gauge invariance of the above Lagrangians and the additional condition that at least one component is uncoloured and electrically neutral, so to provide a viable DM candidate. Considering that the involved SM fields only are quark doublets Q, whose quantum numbers under SU (3) c ⊗ SU (2) L ⊗ U (1) Y are (3, 2, 1/6), and lepton doublets L, (1, 2, −1/2), the new fields have to belong to the representations of SU (3) c and SU (2) L displayed in Table 1. Notice that, while all the possible representations of SU (3) c are listed in the table 1 (as combinations involving larger representations would not feature any colour singlet), in the case of SU (2) L only representations with dimension d ≤ 3 are displayed. The hypercharge assignment is in general not  Table 2 for models with a scalar flavour mediator (Class S).
unique but, as mentioned above, it is restricted by the requirement that at least one state is neutral, i.e. Q = T 3 + Y = 0. Setting the hypercharge of the flavour mediator (Ψ or Φ) as a free parameter X, the hypercharge of other fields are then derived from gauge invariance as shown in the last line of Table 1.
The resulting combinations of quantum numbers are shown in Tables 2 and 3 that collect the possible models with, respectively, a fermion and a scalar flavour mediator. The models have been labelled according to the spin of the flavour mediator (F or S), the combination of SU (2) L and SU (2) c representations of the fields as given in Table 1 and -for categories containing more than one model -the hypercharge of the flavour mediator. The tables include only models featuring at least a viable DM candidate, i.e. an electrically neutral stable state. The representation to which this state belongs have been marked with . Notice that the DM candidate can belong to any of the three NP fields for both classes of models.
DM candidates with non-zero hypercharge are severely constrained by the direct detection experiments, as consequence of the coherently enhanced Spin Independent interactions with nuclei mediated by the Z boson. To keep the particle content of our models minimal, we only consider fermion DM candidates with Y = 0. 2 Instead, if DM belongs to a scalar multiplet one can evade the direct detection bounds by introducing a suitable mass splitting between CP-odd and CP-even components through couplings of the scalar potential. Since this avoids the dangerous DM coupling to the Z without introducing more fields, we include in Tables 2 and 3 also solutions with DM belonging to scalar multiplets with Y = 0. Let us mention 2 Scenarios with Dirac fermion DM and Y = 0 would be still viable if the field content of the model is extended beyond our minimality criterion, such that the DM field mixes with an additional Majorana fermion making the lightest state Majorana, see e.g. Ref. [57,58]. This is the well-known case of the supersymmetric Higgsinos, SU (2)L doublet fermions mixing with a Majorana singlet (the Bino) and a Majorana triplet (the Wino). that some of the models included in the classification shown in Tables 2 and 3 were previously studied in the literature, namely F IA; 0 [34], F IIA [35], S IA [30,33], S IIA; -1/2 [31] and S IIA; 1/2 [32].
In the following sections, we will discuss in details the different constraints (from flavour anomalies, LHC, relic density and direct detection) for the above-defined class of models.

Strategy
In this section we present the strategy that we are going to employ for each model under scrutiny. We take into account constraints coming from flavour physics, from the recasting of direct searches at the LHC and from DM searches. All these constraints are considered in a systematic and comprehensive way in order to assess whether, for a given model, a region of the parameter space where all bounds are evaded exists, or the model is excluded by the combination of all the constraints. We remind that the parameter space of each model is fully determined by three couplings plus three mass parameters, namely for models with a fermion or a scalar flavour mediator, respectively.
Moreover, we also take into account the latest SM results [64,65] In a similar fashion to what we did for b → s transitions, we can interpret the experimental data in terms of a bound to NP effects to a WC: introducing the effective ∆B = 2 Hamiltonian we can observe that the relation between the NP contribution to C BB and ∆M s reads  In order to cast these constraints on the models described in the previous section, we remind that for each of the two classes of models defined in Eqs.
(2)-(3) it is possible to write [28]: and where we have defined Γ Q ≡ Γ Q b Γ Q * s and we have introduced the compact notation The SU (2) L -factors η, η M and η BB are tabulated in Table 4, while the SU (3) c -factor χ M is equal to 1 only if Ψ (Φ) is a Majorana fermion (real scalar) in representation A of Table 1, vanishing otherwise. Moreover, we have introduced the following loop functions The constraints in Eqs. (5) and (8) are exploited through a combined fit to the relevant set of parameters in Eq. (4) using the HEPfit package [66]. For each parameter, we adopt a flat prior with the following ranges: • For the lepton coupling, we allow it to be in the range |Γ L µ | ∈ [0, 4]. Regarding the quark coupling Γ Q we first notice from Eq. (9) that, since the loop functions assume positive values for the mass regimes under scrutiny, it is the only free parameter capable to affect the sign of the Wilson coefficient; hence, in order to obtain the desired sign for δC 9 µ according to Eq. (5) and remembering that Re(V ts ) < 0, we allow it to be in the range Γ Q ∈ [0, 2] for models in class F, and to be in the range Γ Q ∈ [−2, 0] for models in class S; this choice will not affect the result for δC BB , since the coupling is squared in Eq. (10); • For the DM candidate mass, we require for it to be lighter than the other NP fields; • For the remaining NP masses, we let them vary up to 5 TeV. The outcome of the fit is then summarised in terms of posterior probability distribution functions (p.d.f.) for each parameter, together with correlation plots between each pair of them.
As an example, we show in Figure 3 the outcome of the fit to the model F IA; 0 , where the DM candidate is the Majorana fermion Ψ. The diagonal of the triangle plot contains the 1D p.d.f.s for the 5 parameters constrained by the model, while in the other panels we give the combined 2D p.d.f.s showing the correlations among each couple of parameters. Let us now discuss the information we can extract from this result, since it will also be useful for the other models we consider in Section 4. Starting from the couplings 1D p.d.f.s, we observe that Γ Q is constrained to small values: this is a byproduct of the inclusion of ∆M s in the fit, which disfavours higher values for this coupling. However, this implies that |Γ L µ | is pushed to values at the boundary of the perturbative region in order to satisfactory address the B-anomalies. Concerning the masses 1D p.d.f.s, we first notice that the DM mass is allowed to grow only up to a few TeV, due to the requirement that it has to be the lightest state of the NP sector. On the other hand, no upper bound can be inferred on the masses of the other NP fields, which are found to be unconstrained in the whole mass range under scrutiny. Moving on to the analysis of the combined 2D p.d.f.s, we focus our attention on the first column of panels of Figure 3, namely the ones showing the correlations among Γ Q and the other NP parameters included in the fit, where we have highlighted with a dashed red line the value Γ Q ∼ 0.15. From these panels it is indeed possible to infer that such value, sitting in the 1D p.d.f. of Γ Q in the 1σ region and close to its mode, is an interesting benchmark point when looking also at the other parameters. Indeed, from the 2D p.d.f. describing the correlation among Γ Q and |Γ L µ |, we observe that this choice for the quark coupling implies that the 1σ and 2σ regions reach the lowest allowed values for the lepton coupling, corresponding to |Γ L µ | = 3 and 2, respectively, with the 3σ bound reaching |Γ L µ | = 1. Hence, the benchmark value Γ Q ∼ 0.15 is the one that better justifies the benchmark assignments |Γ L µ | = 1, 2, 3, while a different value for the quark coupling would induce bigger values for the leptonic one. Finally, from the other 2D p.d.f.s we observe that in correspondence to this benchmark value the 1σ (2σ) upper bound for M DM is around 700 (1100) GeV, while for the other two fields the whole mass range is allowed at the 2σ level.
An analogous behaviour for the posterior p.d.f.s has been found in all the models analysed in Section 4, with a value for |Γ Q | ∼ 0.15 (positive or negative, depending on whether the model belongs to class F or S respectively) such that the allowed 1σ and 2σ bounds on |Γ L µ | are the lowest, the 1σ and 2σ upper bounds for the DM mass are around 700 GeV and 1100 GeV respectively, and no upper bound is found at the 2σ level on the remaining masses. Therefore, in the following analysis we will fix the value for the quark coupling to the benchmark value |Γ Q | = 0.15, and we will consider 3 different benchmark values for the leptonic coupling |Γ L µ | = 1, 2, 3, corresponding to the 3σ, 2σ and 1σ lower bounds, respectively. Moreover, we observe that the fit always allows to fix the mass of one the two heavy NP partners of the DM candidate to be equal to 700 GeV (1100 GeV), value corresponding to the 1σ (2σ) upper bound for the DM (which being the lightest NP particle obviously sets a lower bound on the rest of the spectrum).
Let us finish this subsection with a comment on another observable featuring a long-standing deviation from the SM prediction: the muon anomalous magnetic moment, (g − 2) µ . Notice that the subset of our NP fields coupling to muons indeed contributes to the dipole operator responsible for (g − 2) µ . Such a contribution is, however, chirally suppressed (our fields do not couple to right-handed muons), thus too small to account for the current ∼ 3.5 σ discrepancy between SM prediction and experimental measurement [67], with the possible exception of tuned regions of the paramater space, as systematically discussed in [68]. A successful fit of the (g − 2) µ discrepancy would require the introduction of extra fields coupling to the SM Higgs field and thus providing a chiral enhancement to the muon dipole. Examples of such "next-tominimal" dark matter scenarios naturally accounting for the (g − 2) µ discrepancy are given in [68,69] as well as, in combination with the B-physics anomalies, in [50,51].

LHC constraints
The second step of our strategy involves studying the bounds following from direct searches at the LHC. We will consider here as an example, without any loss of generality, the case where the flavour mediator is a scalar field (being the conclusions in the case of a fermion mediator identical). DM scenarios and LHC searches can be separated in two categories, depending on whether the DM candidate is the fermion field coupling to SM leptons (Ψ L ), or one of the two other possibilities (Ψ Q or Φ).
• Case 1. If the DM candidate is either the field Ψ Q or Φ, the main production channel at the LHC will be via QCD-mediated processes, such as with DM being the field Φ or Ψ Q , respectively, and appearing as missing energy in the detector. Here q, q = s, c, b, t, that is, each of the quarks can be either a light quark, producing a jet, a bottom quark, producing a b-jet, or a top quark, producing a b-jet plus the products of the decay of the W boson. 3 Notice the similarity of this setup to simplified supersymmetric (SUSY) models with squarks decaying into neutralino DM. As a consequence, in order to produce exclusion plots in the (M Ψ Q , M Φ ) plane, it is possible to recast limits from LHC squark searches involving missing energy plus 2-6 jets, allowing also the presence of b-jets. The latest of such analyses has been performed by CMS [70]. Let us stress once again that this recasting can be applied no matter which of the two NP fields in the couple {Ψ Q , Φ} is the DM candidate, since the final states and hence the LHC signature is the same for both cases, as shown in Eq. (12).
In order to set bounds on M Ψ L , we now have to distinguish whether the DM candidate is the scalar Φ or the fermion Ψ Q coupling to SM quarks. In the former case, one can constrain M Ψ L with a procedure analogous to the one outlined above: a pair of Ψ L fields will be mainly produced by electroweak Drell-Yann, and they will subsequently decay into muons and DM, i.e. missing energy. Hence, this channel can be described as Therefore, similarly to what described above, in order to produce exclusion plots in the plane (M Ψ L , M Φ ) it is possible to recast limits from LHC slepton searches involving missing energy plus a muon pair. The latest of such analyses has been performed by ATLAS [71]. Moreover, there are also searches considering soft leptons performed by both ATLAS [72] and CMS [73]. In the second case (DM in Ψ Q ), on the other hand, Ψ L does not couple directly to the DM candidate, hence the above signature cannot be used. The main decay channel of the Ψ L fields will consist instead in the cascade production of two muons plus two NP scalars Φ, which will further decay into quarks and DM candidates, i.e. missing energy. Hence, in this last scenario the signature will be Therefore, in order to produce exclusion plots in the plane (M Ψ L , M Ψ Q ), it is necessary to recast limits from combined LHC searches for stops cascade decaying into sleptons/charginos, involving missing energy plus 2 muons and 2-6 jets (including b-jets). The latest of such analyses can be found in [74,75].
• Case 2. If the DM candidate is the field Ψ L , analogous considerations to the ones reported above Eq. (13), can be applied, with the signature reading now One can therefore use again the searches from Refs. [71][72][73] to constrain the plane (M Φ , M Ψ L ). Further, bounds on M Ψ Q can be obtained through the cascade decay Hence, exclusion plots in the plane (M Ψ Q , M Ψ L ) can be obtained again by means of the analyses in [74,75].
While some combination of the above production and decay modes will appear in all models under study, there is a further collider signature that is possible if DM is part of an SU (2) L multiplet, where charged states will be also present. Due to electroweak radiative corrections the charged states are typically O(100) MeV heavier than the neutral DM state [57,76,77]. Because of such a small mass splitting, the charged states -that can be pair produced through electroweak Drell-Yann at the LHC -are long-lived, that is, they can travel a macroscopic distance (typically a few cm) in the detector before decaying (through the exchange of an off-shell W boson) into DM and a very soft and thus undetectable pion. The resulting signature is a so-called "disappearing track" observable in the inner tracker of the detector. Searches for this kind of events have been performed by both ATLAS [78] and CMS [79].
In order to exploit the LHC searches discussed above 4 we have implemented, as first step, the models in Feynrules 2.3 [80] in order to generate UFO files, which are subsequently passed to MadGraph5 aMC@NLO 2.7 [81] where the matrix elements are calculated and a set of 50k events is generated. The partonic events are showered using Pythia 8 [82], the detector effects are simulated by means of DELPHES 3 [83] and the result is eventually passed to CheckMATE 2 [84], which compares the simulated signal with the experimental searches at the LHC and determines whether the model point is excluded at the 90% confidence level. For the LHC searches selected above that are not yet implemented in the current version of CheckMATE, we employ the AnalysisManager framework in order to define them ourselves [85]. As said above, such procedure determines whether a model is excluded at the 90% confidence level for a fixed value of the 6 This means that, in order to produce an exclusion plot in the (M Ψ Q , M Φ ) plane (or in the (M Ψ L , M Φ ) one, if Ψ L is the DM candidate), we have to fix the values of the 3 couplings and of the remaining mass to some set values. The chosen values for these 4 parameters will be guided by the benchmark values inferred from the flavour fit performed in the previous section. Regarding the lepton coupling, we can directly adopt one of the 3 benchmark values inferred from the flavour fit, i.e. |Γ L µ | = 1, 2, 3. However, the same cannot be done for the quark couplings since the flavour fit gave a benchmark value for the product |Γ Q b Γ Q * s |, i.e. |Γ Q | = 0.15. On the other hand, LHC constraints are sensitive to the individual values to the two quark couplings and not only to their product. Hence, in the following we will inspect for each model 2 different benchmark cases: We are now left to fix the value of one of the 3 masses. In the case of the DM being Ψ L , we will fix the value of the mass of Ψ Q , which is a coloured particle not directly coupling to the DM candidate. The produced signature will be the cascade decay shown in Eq. (16), and hence we will have to rely on the recasting from Refs. [74,75]. We infer that (unless M Ψ Q − M Ψ L 50 GeV) we have to consider as a benchmark value M Ψ Q = 1400 GeV. This is indeed the lowest allowed value for this parameter such that no lower bound is induced on the DM mass. In a similar fashion, if the DM candidate is Ψ Q we will fix the mass of the coloured field Ψ L . Once again, given the signature of Eq. (14) we employ the recasting from Refs. [74,75] and obtain also for this field the benchmark value M Ψ L = 1400 GeV. Finally, if the flavour mediator Φ is the DM particle, the particle whose mass we will fix is yet again Ψ L , but this field is now a colour singlet directly coupling to the DM candidate. We can therefore directly employ the recasting from Refs. [71][72][73], and observe that we can set the mass of Ψ L to the one of the 2 benchmark values inferred from the flavour fit of the previous section, i.e. M Ψ L =700 GeV or 1100 GeV. Indeed, the highest excluded value from the recasting of such results is always found to be at most O(600) GeV.
Summarizing, given the above benchmark values for 4 parameters, it is possible to combine either in the (M Ψ Q , M Φ ) plane or in the (M Ψ L , M Φ ) one (according to whether the DM candidate is Ψ Q or Φ, in the former case, or Ψ L , in the latter). Indeed, the WC δC 9 µ = −δC 10 µ is simply a function of the 2 masses, once all the other parameters are fixed to benchmark values. Following the procedure outlined above it is therefore possible to visualise on the same plane both the region that can account for the flavour anomalies at the 2σ level as in Eq. (5), and the one excluded at 90% CL by our recasting of direct searches at LHC.

Constraints from DM phenomenology
As a last step, the flavour and LHC constraints can be complemented with the ones from DM phenomenology. Concerning the latter, the models considered here belong to the category of the so-called t-channel portals [86][87][88][89][90][91] (a "flavoured" variant of this kind of setup has been also considered here [92][93][94][95][96]). Our assumption that a good fit of the B-anomalies is achieved introduces, however, some relevant variation in the phenomenology of this kind of models, especially for what concerns direct detection.
As will be shown in the following, being some model parameters fixed by the fit of flavour observables, the latter constraints can be easily visualised into two-dimensional mass-mass plots to compare the corresponding viable regions with those fulfilling the requirements from both LHC and flavour physics.

Relic density
The DM relic abundance has been measured with great precision by the Planck experiment [97] and it is represented by the parameter Ω DM h 2 whose value is: As will be evident in the next sections, for values of the couplings compatible with flavour anomalies, the DM is capable of reaching thermal equilibrium in the Early Universe and, hence, can achieve its final relic density through the freeze-out mechanism. In such a case, the DM abundance is the solution of a Boltzmann equation which can be written as [98,99]: where σv eff is the effective thermally-averaged DM annihilation cross-section including coannihilations and g * is the effective number of relativistic degrees of freedom. The integral is computed between the freeze-out temperature ) and the current temperature of the universe T 0 . Defining M the field that is t-channel exchanged in DM pair annihilations and also contributes to co-annihilation processes, the DM effective annihilation cross-section can be written as [86]: where x = M DM /T . σv DM DM describes DM pair annihilation processes into SM fermions mediated by t-channel exchange of the field M . Given the assumptions stated in the previous sections, the possible final states are µ + µ − , andqq where q, q = s, c, b, t. 5 σv DM M represents coannihilation processes with a DM and a M particle in the initial states while σv M † M describes the contribution of M pair annihilation processes to the DM effective annihilation cross-section provided that the mass splitting between DM and M is sufficiently small. Notice that the expression above is valid for complex scalar and Dirac fermion DM.
In the case of real scalar and Majorana DM we have a slightly different expression: In the above equations ∆ = (M M − M DM )/M DM is the relative DM/mediator mass splitting while: where g M and g DM are the internal degrees of freedom of the mediator and of the DM. In our analysis, the DM relic density, including coannihilations, have been computed with great numerical precision through the package micrOMEGAs [37]. To clarify our results we provide nevertheless analytical expressions of the DM annihilation cross-section into fermion pairs, the most relevant in the regions of parameter space favored by B-anomalies (see below), at the leading order in the conventional velocity expansion (as further simplification we have neglected the masses of the final state fermions) [86][87][88]100]: where N c = 3(1) in case of colour charged (colour singlet) final state fermions. The sums run over the kinematically accessible final states, depending on the value of the DM mass.
The four expressions refer, as indicated, to real scalar, complex scalar, Dirac fermion and Majorana fermion DM. Scalar and fermionic DM candidates have been generically called Φ DM and Ψ DM , respectively, while the states exchanged in the t-channel Feynman diagrams and interacting with the fermion f have been called S f and F f . Finally λ f correspond to suitable assignments of Γ L,Q µ,s,b according to the final states. 6 While the four cross-sections have a very similar mass dependence in the limit in which DM is much lighter than the NP field exchanged in the t-channel, they feature a very different velocity dependence. We notice indeed that the annihilation cross-section is s-wave dominated, in the case of Dirac fermion, p-wave suppressed in the case of complex scalar DM and Majorana fermion, and even further (d-wave) suppressed in the case of real scalar DM. Given that v 2 ∼ 0.1, we expect that, while Dirac DM will easily comply with the 5 Up quarks can be as well annihilation final states. Their contribution to the DM relic density is, however, negligible because of the CKM suppression of the couplings. 6 Notice that colour charged NP field might interact with different quark generations. The expression of the annihilation cross-section might slightly change in this case.
requirement of the correct DM relic density, real scalar DM will, instead, by typically overabundant in light of its very suppressed annihilation cross-section, unless the latter will be enhanced e.g. by coannihilations.
Annihilations into SM fermion pairs and coannihilations mediated by M represent the main contribution to the DM relic density in the case the DM belongs to an SU (2) L singlet. In the case DM belongs to an SU (2) L multiplet, it can also annihilate, through gauge interactions, into W and Z boson pairs. The latter annihilation processes easily become the dominant contribution to DM pair annihilations since the corresponding annihilation rate is not suppressed by the mass of the field M . For DM masses above the TeV, such cross-section is further increased by the so called Sommerfeld enhancement [101][102][103][104] as well as by bound state formation [105]. Additional coannihilation processes, due to other components of the DM multiplet, are present as well.
As will be clear from the following analysis, imposing the correct relic density, Eq. (18), translates into a very strong constraint, only marginally compatible with the flavour anomalies and other phenomenological bounds. For this reason we will just apply, through Eq. (19), an overclosure bound Ω DM h 2 ≤ 0.12. However, while requiring that thermal DM production does not exceed the observed relic density, we will also assume that our DM candidate always accounts for 100% of DM, including in the regions of the parameter space where it would be underabundant, as a consequence of some (unspecified) non-thermal DM production mechanism, see e.g. [106,107].

Direct detection
The scattering of the DM with nucleons and nuclei, which is at the base of direct detection (DD), is typically described through effective four-field operators coupling pairs of DM particles with SM quark or gluon pairs. For all the models considered in this work, the strongest constraints come from Spin Independent (SI) interactions. For our analysis we have adopted the world leading limits given by the XENON1T collaboration [108]. How effective are the resulting constraints depends on the spin of the DM and, in the case of scalar DM, on whether the field is real or complex, while in the case of fermionic DM, on its Dirac rather than Majorana nature. In the following illustrative discussion, we focus for simplicity on the case in which the DM is an SU (2) L singlet.
In the case of complex scalar DM, the effective Lagrangian for DD reads: where O q µν and O g µν are the twist-2 operators: As we can see, the effective Lagrangian considers just interactions with light quarks (q = u, d, s) and gluons. This is because the typical energy scale for DM scattering processes is of the order of 1 GeV and, hence, heavy quark flavours, c, b, t, are integrated out.
The coefficient c q in Eq. (23) can be further decomposed as c q = c q tree + c q Z + c q γ + c q box . Illustrative diagrams associated to these different contributions are shown in Figure 4. Notice that the figure actually displays the case of fermionic DM, discussed in the following. However topologically analogous diagrams can be also obtained for scalar DM, since as we will discuss below, the two cases share many similarities. c q tree is the tree-level induced contribution from diagrams with s-channel exchange of colour charged scalar NP field. Since such tree level contribution describes the interactions of vector currents, the corresponding coefficient for the nucleons is just a linear combination of the contributions of the valence quarks, namely: In the models considered in this work the c N =p,n tree coefficients are generated only be the CKM mixing. Having chosen a basis in which the down-type quark mass matrix is flavour diagonal we have: Given this result, one cannot neglect a priori contributions from loop-level induced interactions. Indeed, the c q Z,γ coefficients are generated by penguin-like diagrams as the ones shown in the second panel of Figure 4, charged under the SM EW group, with SM γ and Z bosons while c q box is the coefficient associated to box diagrams analogous to the one present in the third panel of Figure 4. The remaining operators in Eq. (23) arise, again at the loop level, from QCD interactions of the colour charged new fermions, possibly present in the theory [109]. Full analytical expressions for the c q Z,γ , c q box , d q,g , g q,g 1 can be found e.g. in Refs. [30,109,110]. As will be discussed in the following, the strongest limits will arise from γ, Z penguins. We report here the relevant simplified expressions, as given in [30], in the limit For what concerns the γ penguin we have: The coefficient associated to the Z-penguin, in the same approximation, instead reads: 7 where Among these two contributions, photon penguins give typically the dominant contribution with the exception of the case in which the effective coupling of the DM with top quarks is sizable, as a consequence of the enhancement proportional to the square of the top mass in c q Z . The effective Lagrangian in Eq. (23) gives rise to the following scattering cross-section for the DM over nucleons (for illustration we focus on the case of the proton): where µ p = M Φ DM M p /(M Φ DM + M p ) is the DM/proton reduced mass. The extra factor depending on A, Z, being respectively the mass and atomic number of the detector material, allows a consistent comparison with experimental limits which assume equal coupling of the DM with protons and neutrons [111]. f p,n represent, in fact, the effective coupling of the DM with protons and neutrons and read: where c p i = 2c u i + c d i , c n i = c u i + 2c d i with i = tree, Z, γ, box, and f G = d g c + d g b + d g t . The parameters f N =n,p q , f T G , q(2), and G(2) are nucleon form factors defined as: For our analysis we have used the default values implemented into the micrOMEGAs package [37].
It is important to remark that, contrary to the other coefficients, including the ones generated at the tree-level, c p γ,Z coefficients can be present in models in which the DM is coupled only to NP states charged under the electroweak gauge groups but that are colour singlets.
If DM is a real scalar, the Φ † DM i ↔ ∂ µ Φ DM operator identically vanishes. Hence, the DM direct detection cross-section is expected to be suppressed.
Moving to fermionic DM, which we just call Ψ DM , the effective Lagrangian for the Dirac case is: Again, the coefficient c q is a combination of CKM-suppressed tree-level and loop-induced contributions: where, in contrast to the case of scalar DM, a contribution from Higgs penguin diagrams is present as well. The operatorΨ DM γ µ Ψ DMq γ µ q behaves, with respect to direct detection, in an analogous way as So we have again that the coefficients c p,n at the nucleon level are linear combinations of the coefficients associated to up and down quarks. The coefficients d q,g can be decomposed into QCD contributions, which we label d q,g QCD , analogous to the ones discussed for scalar DM, and a contribution from Higgs penguin diagrams, labelled as d q,g H . The coefficients g g,q 1 finally come from QCD interactions mediated by quarks/gluons and possible NP colour-charged states. Effective interactions mediated by the photon are present as well. The latter are described, this time, by the following Lagrangian: with the two terms dubbed, respectively, magnetic dipole moment and charge radius operators.
Concerning the relative contribution of the different coefficients, as already discussed in Ref. [112], the situation is analogous to the case of complex scalar DM. The dominant contribution is typically associated to the charged radius and dipole operators, whose coefficient can be approximately written as: Again, in case of sizable couplings of the DM with the top quark, the latter terms are overcome by c q Z which can be approximately written as: Given the presence of dipole operators, direct detection phenomenology is not fully caught by the scattering cross-section over nucleons but, on the contrary, one has to rely on the DM scattering rate over nuclei: where and µ T = M Ψ DM M T /(M Ψ DM + M T ) with M T being the mass of target nucleus. F SI is the conventional SI nuclear form factor [113], while F D is the form factor associated to dipole scattering [114]. Experimental limits have been obtained, in this case, with the procedure illustrated in Ref. [115].
In the case of Dirac DM the coefficients f p and f n are written as: Changing the nature of the DM, this time from Dirac to Majorana, leads to markedly different case. Indeed, theΨγ µ Ψ andΨσ µν Ψ operators are identically null and the DD phenomenology is again fully captured by the SI scattering cross-section of the DM on protons: where f p,n are defined analogously to Eq. (40). In the case in which the DM belongs to an SU (2) L multiplet, no new operators are generated in the Lagrangians in Eqs. (23) and (33). The coefficients of the effective operators get, however, additional contributions from loop diagrams in which Z, W bosons are exchanged. The case of Majorana DM has been discussed extensively e.g. in Refs. [58,116] while the case of real DM has been considered, to a lower extent e.g. in Ref. [109]. To our knowledge no analogous computations are available for complex scalar and Dirac fermionic DM.

Indirect detection
As well known, indirect detection (ID) for WIMPs, relies on the search of the products of residual annihilation processes for DM occurring at present times. Similarly to the case of direct detection it is convenient to distinguish the cases in which the DM is an SU (2) L singlet or not. In the former case, the main annihilation channels to consider are the ones into SM fermion pair final states. These lead mostly to continuous γ-ray signals which, for the ranges of DM masses considered in this work, can be probed by telescopes such as Fermi-LAT [117,118]. The impact of the resulting constraints is highly model-dependent though. Indeed, having in mind the velocity expansion: σv ≈ a + bv 2 , we have that only for s-wave dominated annihilation cross-section, i.e. a = 0, the values of the cross-section at thermal freeze-out and and present times are comparable, so that eventual ID limits are effective. On the contrary, p-wave, i.e. the b coefficient is the leading contribution, dominated cross-section are affected by ID to a negligible extent. Notice as well that coannihilations are also mostly effective at thermal freeze-out while their rates are, instead, exponentially suppressed at present times. Given this, among the models presented in this work, only scenarios with Dirac fermionic DM can be probed by indirect detection. 8 Summarising, in the following section we will apply the strategy here described to several models of interest, in order to study whether such models allow for a region of the parameter space where all the constraints here described are evaded, or the model is excluded by the combination of all the constraints.

Results and discussion
In this section, we analyse some of the models listed in Tables 2 and 3, following the strategy illustrated in the previous section. Such models have been highlighted in cyan in the tables. Our selection covers a broad variety of cases, including both scalar and fermion flavour mediators, as well as both scalar and fermion DM. Furthermore we will separately discuss, where appropriate, both real and complex scalar DM as well as both Dirac and Majorana nature for fermionic DM. Finally, notice that our selection comprises models where DM is a pure SM singlets as well as cases where it belongs to SU (2) L multiplets.
We will write for each model the Lagrangian responsible for the phenomenology we are interested in, according to the quantum numbers of the NP particles, and determine the regions of parameter space for which the B-anomalies are accounted for. We will then combine this requirement with the constraints from collider searches of the NP particles as well as from DM phenomenology, in particular relic density and direct detection. As discussed in Section 3.3.1, we will assume that in the regions of the parameter space where thermal DM production is insufficient some non-thermal mechanism is at work such that our DM candidate always account for 100% of the observed DM abundance.

F IA; 0 , Dirac singlet DM
We start considering the model F IA; 0 with singlet Dirac DM. This case, which is among the simplest in Tables 2 and 3, has never been studied before and, as we will see, is subject to strong constraints. It is a good example to illustrate how bounds from different sources can altogether exclude a model. The Lagrangian of this model reads: with the fields Φ Q , Φ L and Ψ carrying respectively the following SU (3) c ⊗ SU (2) L ⊗ U (1) Y quantum numbers: (3, 2, 1/6), (1, 2, −1/2) and (1, 1, 0). As mentioned above, the DM candidate is the Dirac field Ψ DM = Ψ, which also plays the role of the flavour mediator in the diagram in Fig. 1. This scenario has been selected for our analysis since it features the highest degree of correlation. Indeed, notice that the DM field couples to both the NP fields Φ Q and Φ L , which are charged under the SM gauge group. Consequently, all the three couplings Γ Q s , Γ Q b , Γ L µ , entering δC 9,10 µ are, as well, contributing to the DM annihilation and scattering rates. For this reason we will also investigate, in the following, similar models with Majorana fermion and scalar DM.
The results of our analysis for this model are presented in Figure 5, displayed in the (M Ψ , M Φ Q ) twodimensional plane. The two rows in the figure correspond to the two different assignments of the couplings , which proved to provide an equally good fit of the B-anomalies, cf. Section 3.1. For each coupling configuration we show only one of our benchmark values for M Φ L , namely 1100 GeV, and three values of |Γ L µ |. In each plot, the region compatible with the flavour anomalies is the one enclosed within the two green contours. As can be seen from the different filling styles, the regions outside these bands should be interpreted differently. Indeed the regions on the left of the green contours (filled in green) are ruled out, since they correspond to the case in which δC 9,10 µ exceed the experimental limits. On the right of the contours, on the contrary, NP contributions to C 9,10 µ are increasingly suppressed so that the these observables do not deviate, to a statistically relevant extent, with respect to the SM expectation. While current flavour anomalies are not reproduced in the latter parameter regions, we cannot strictly regard them as ruled out as the anomalies are still awaiting full experimental confirmation. These regions are denoted by a green horizontal hatching.
The orange region represents the exclusion from LHC searches for the signatures with jets and/or muons and missing energy described in Section 3.2. For this model, we show our recasting of the bound from Ref. [70] Moving to DM phenomenology, the constraints from DM relic density are represented as red regions. As already mentioned, throughout our study we will just require that the value of Ω DM h 2 , determined by applying the conventional thermal freeze-out paradigm, does not exceed the experimental determination, namely Ω DM h 2 ≤ 0.12. In each plot the region of parameter space which does not fulfill this constraint has been marked in red. Being the DM a SM singlet, its relic density is determined, with the exception of the coannihilation region, by annihilations into SM fermion pairs. As the value of Γ L µ increases, the region with overabundant DM progressively reduces and, for Γ L µ = 3 the DM is always underabundant within the whole range of M Ψ and M Φ Q shown in the plot. The blue region corresponds, instead, to the case in which the DM interactions with nuclei, as given by Eq. (38), exceed the constraints from XENON1T. Finally, being the DM a Dirac fermion, constraints from indirect detection should be taken also into account. The regions of parameter space at the left side of the dashed yellow contour are excluded by the latter type of searches.
As evident, in all the plots the region compatible with the flavour anomalies falls at least into one of the experimental exclusions. Among them, the strongest by far comes from direct detection, which excludes the whole range of masses considered in the different plots, besides the case Γ L µ = 1: with such an assignment, only a region with M Φ Q = 4 − 5 TeV survives in case i ), while a broader area is allowed in case ii ). Nevertheless, in these the NP contribution δC 9 µ = −δC 10 µ is too small to account for the observed anomalies. The direct detection bound extends even beyond multi-TeV masses for the Φ Q field because it is actually saturated by the charge radius and magnetic dipole operators in Eq. (38) which are dominated by the contribution of the colour singlet field M Φ L , whose mass and coupling Γ L µ were kept fixed in the analysis. The case M Φ L = 1.1 TeV is hence ruled out by DM DD regardless the assignment of the other parameters. For the same reasons we have shown no plot for M Φ L = 700 GeV since, in such a case, also the case Γ L µ = 1 would be completely ruled out by direct detection.

F IA; 0 , Majorana singlet DM
In this subsection we study the same model discussed in the previous one, defined by the Lagrangian in Eq. (43). The only difference is the nature of the DM field Ψ, now corresponding to a Majorana fermion. As pointed out above, this kind of scenario is particularly interesting since it features the highest degree of correlations between B-anomalies and the other phenomenological observables considered in the present study. A model analogous to F IA; 0 with Majorana DM has been already studied in Ref. [34]. The analysis in the latter reference differs from the present work in the fact that they strictly imposed the requirement of the correct relic density and used it to fix M Φ L as a function of the other parameters, in particular the coupling Γ L µ . Furthermore, the latter parameter has been allowed to reach the perturbativity bound Γ L µ = √ 4π. As already mentioned, we adopt here a more conservative approach for what concerns the DM relic density, keeping M Φ L and Γ L µ as free parameters in the fit of flavour observables and considering Γ L µ ≤ 3 in our phenomenological study. Given these different assumptions our findings slightly differ from the ones reported in Ref. [34]. We repeat the analysis whose procedure has been illustrated in detail in the previous section. The results are shown in Figures 6 and 7. As we can see, we find a very different picture with respect to the results of the previous case as shown in Figure 5. This is a consequence of the different nature of the fermionic DM candidate. One can indeed find regions of the parameter space compatible with the flavour anomalies and not in tension with other experimental bounds, even for the lighter benchmark value M φ L = 700 GeV. For Γ L µ 2 it is also possible to have, albeit in a narrow window of the parameter space, a good fit of the flavour anomalies and, simultaneously, saturate the observed DM relic density (Ω DM h 2 0.12) with a standard thermal WIMP. This occurs for DM masses between approximately 50 and 150 GeV.
This different outcome compared to the Dirac case is mostly due to the fact that, for Majorana DM, the effective operators accounting for DM interactions with nucleons mediated by the Z boson and the photon identically vanish. As a consequence, the strength of the interactions between DM and nuclei is strongly reduced and, hence, DD bounds are significantly weaker. Furthermore, the annihilation cross-section of Majorana DM is p-wave suppressed, and thus ID constraints are not present. Although this is a successful scenario according to our criteria, we remark again that, in order to have compatibility between the fit

F IB; -1/3 , Real-scalar doublet DM
We now turn to consider an example with scalar DM. In this case, DM is a (real) neutral state which is part of the scalar SU (2) doublet Φ Q , (1, 2, 1/2). The other NP fields Φ L and Ψ transform under the SM gauge group respectively as (3, 2, −1/6) and (3, 1, −1/3). The Lagrangian resembles that of the previously-considered model: This is the simplest model, with our general classification, featuring the DM belonging to a SU (2) doublet. As discussed in Section 2, our minimality criteria only allow scalar DM to belong to an SU (2) multiplet. Notice also that, as usual, the DM field Φ Q and the other NP scalar Φ L could couple with the Higgs field and among each other through operators of the type |Φ Q | 2 |Φ L | 2 , |H| 2 |Φ Q | 2 and |H| 2 |Φ L | 2 . As already pointed out, we are assuming in this work that these quartic couplings are so small that have a negligible impact on phenomenology. The only exception is given by the coupling providing a mass splitting between real and imaginary parts of the neutral component of Φ Q , as discussed below. In fact, since the DM belongs to an SU (2) doublet with non-zero hypercharge, it would feature as well tree level interactions with the Z-boson. In order to circumvent this possibility (which is already experimentally ruled out, see e.g. Ref. [88]), we assume that, similarly to the so-called inert doublet model [119][120][121][122], the neutral component of Φ Q can be separated into a CP-even state, which we assume here to be the DM candidate, and a CP-odd state with a sufficient mass splitting ( O(100) keV) to avoid DD constraints. This can be achieved through a quartic operator involving the SM Higgs of the form (|Φ † Q H| 2 + h.c.). Since a tiny coupling O(10 −13 ) is enough to induce an O(100) keV mass splitting, we can safely assume that no other phenomenologically relevant effect follows from this (and other) Higgs-portal interactions.
The combined constraints on this model are shown, with the usual colour coding, in Figures 8 in the (M Ψ , M Φ Q ) plane (notice that what we label M Φ Q is just the DM mass). In contrast to the previous model, the strongest constraints come from flavour and LHC. 9 First, notice that we set M Φ L = 1400 GeV, in order 9 Here we have not considered bounds from disappearing tracks that would arguably have little sensitivity to the small to evade bounds on Φ L production and cascade decay (pp → Φ L Φ L → µ + µ − + ΦΦ → µ + µ − + qq + / E T ) from the LHC searches in Refs. [74,75]. The dominant collider bound shown in the plots follows from the process pp → ΨΨ → qq + Φ Q Φ Q , to which the CMS search of Ref. [70] is sensitive. On the other hand, the flavour anomalies are accounted for only for relatively light values of M Ψ which fall in the region excluded by the CMS search, 10 with the exception of narrow tuned strips on the (M Ψ , M Φ Q ) plane corresponding to a very compressed spectrum, where jets would be too soft for a substantial number of events to be selected by the experimental cuts. We remark again that, even if it is not capable of accounting for the B-anomalies compatibly with all experimental constraints, the model F IB; -1/3 is not strictly ruled out since LHC and DM constraints are evaded in large regions of the parameter space where the NP contributions to B-observables does not exceed, in a statistically relevant way, the SM ones. Moving to DM phenomenology, we can see from Figure 8 that direct detection constrains the model under consideration only poorly. This occurs because of the nature of DM: besides excluding tree-level interactions between DM pairs and the Z boson, the small mass splitting we assumed also sets to zero the Φ † Q i ↔ ∂ µ Φ Q operator in the DD effective Lagrangian in Eq. (23). Finally, concerning DM relic density, we note that, similarly to what occurs in the case of the inert doublet model [119][120][121][122], very efficient annihilation processes into gauge bosons dominate DM production, as long as they are kinematically accessible. Thus relatively large values of the DM mass (≈ 500 − 600 GeV) outside the region fitting the B-anomaly are required to have the correct relic density.

S IA , Complex-scalar singlet DM
This model is the counterpart of F IA; 0 with scalar and fermion fields exchanging roles. It is then of primary interest due to the high degree of correlation among our observables. While a similar study of this scenario has been presented already in Ref. [30], our results are notably different, as discussed below.
The DM candidate belongs to the complex field Φ, which is a complete singlet under the SM gauge group. Since Φ also plays the role of the flavour mediator in the diagram of Fig. 1, it is not possible in this case to assume that it is a real scalar field, otherwise an additional "crossed" box diagram would exactly cancel the effect on b → s , as can be seen from Eq. (9). The quantum numbers of the other two fields, Ψ Q and Ψ L , are respectively (3, 2, 1/6), (1, 2, −1/2), and the Lagrangian of the model reads: Given our results for the F IA; 0 model discussed in Section 4.1, since the Φ † i ↔ ∂ µ Φ operator behaves, in the non-relativistic limit relevant for DM DD, in an analogous way as the Dirac DM operatorΨγ µ Ψ, we can expect this model to be completely ruled out by DM direct detection as well. Furthermore, as we mentioned above, we can not assume Φ to be a real scalar if we want the model to address the B-physics anomalies. To circumvent this problem, we can however assume that the two components of Φ (the real part and the imaginary part) have a small mass splitting > O(100) keV. 11 In this "non-degenerate" case the DM candidate would be a real scalar, the lightest of the two states. While this would ensure the needed suppression for the scattering cross-section with nuclei, a small mass splitting would avoid the production cross section and life-time of the charged states of the doublet. 10 This bound extends to substantially larger masses than in the previous case FIA; 0, as the produced particle here is a fermion, hence its production cross section is about one order of magnitude larger than that of a scalar of the same mass. 11 If DM is stabilised by a Z2 symmetry, this can be achieved for instance through an invariant term in the scalar potential of the form µ 2 (Φ 2 + Φ * 2 ) [35].  Figure 5. In particular, the regions strictly excluded by DD are filled in blue, while the areas denoted by a vertical blue hatching are excluded only if the two components of the DM singlet are degenerate, see the text for details.
above-mentioned cancellation of δC 9,10 µ (the DM is the flavour mediator in this case) and maintain efficient DM annihilations. In fact, the annihilation cross-section into SM fermions would be d-wave suppressed for real DM, hence the correct relic density would follow from coannihilations of the two components of Φ.
The scenario with two degenerate states has been studied in Ref. [30]. 12 Considering also the nondegenerate case will allow us to open a large viable region of the parameter space. Besides this important point, Ref. [30] differs from our study also from the fact that LHC data available at the time lead to less stringent bounds for couplings as in scenario ii ) compared to the limits in scenario i ). This is however not the case anymore, since the two scenarios are now no longer significantly different once the recasting of the CMS search of Ref. [70] is employed.
The combined constraints for this model are shown in Figures 9 and 10. In these plots, the DD constraint 12 A slightly more complicated variant of this model, featuring an additional scalar field, has been studied in Ref. [33].
Furthermore, adding to the model a scalar doublet mixing with the singlet through a coupling with the Higgs introduces interactions to right-handed fermions as well, allowing in particular a chirally-enhanced contribution to (g − 2)µ and thus a natural explanation of the observed anomaly [50].

S IIB , Dirac singlet DM
We now discuss another model with fermion singlet DM. The gauge quantum numbers of the NP fields are (1, 1, 0), (3, 1, −2/3) and (3, 2, 1/6) for Ψ Q , Ψ L and Φ, respectively. The DM candidate is the Dirac field Ψ Q , and the flavour mediator Φ is a complex scalar: This model is characterised by the fact that the DM field only couples to quarks and not to muons. It then allows for an interesting comparison with the F IA; 0 model. Similarly to the latter scenario we will consider, respectively here and in the next subsection, both cases of Dirac and Majorana DM. The results of our analysis are shown in Figure 11 and confirm the tendency, already observed for F IA; 0 , of models with Dirac DM of being already experimentally excluded. Compared to the F IA; 0 scenario we notice some differences though. First of all the relic density constraint is much stronger and, moreover, is not sensitive to the value of the Γ L µ coupling. In fact, all DM phenomenology does not depend on this latter parameter (neither on M Ψ L ), as it can be seen by comparing the three columns of plots of Figure 11. This is due to the fact that the DM field is now only coupled with the field Φ and left-handed quarks. DM observables hence depend only on the Γ Q s,b couplings which are fixed by the fit of flavour observables. Furthermore, being the coupling with muons absent, DM annihilation cross-section is more suppressed than in the F IA; 0 model, especially for the configuration ii ) of the couplings. Looking instead at direct detection, while being still the most challenging constraint for the assignment i ) of the couplings, its impact is strongly reduced, relative to the requirement of viable relic density, for the assignment ii ). This is because the most relevant contribution to the SI cross-section here come from loop diagrams involving the Z bosons, whose effect is enhanced by the mass of the top quark. Since the coupling of the DM field with the top is reduced, once moving from the assignment i ) to ii ), the DM scattering rate is consequently reduced. We have not observed such an outcome in the model F IA; 0 since a comparable or even larger contribution to the DM scattering rate on nucleons was coming from interactions with the photon, controlled by the Γ L µ coupling. As mentioned above, in the S IIB models, all the phenomenological constraints are independent on the M Ψ L and Γ L µ parameters (as long as Ψ L is heavy enough to evade searches for cascade decays as in Eq. (14)) with the exception of the region favoured by B-physics. Moreover we have found that the latter changes only marginally for different values of M Ψ L . Hence we have set for all plots M Ψ L = 1400 GeV, a value that allows to evade the constraint on pp → Ψ L Ψ L → µ + µ − + ΦΦ → µ + µ − + qq + / E T of the searches in Refs. [74,75]. Figure 11 shows that, mainly due to the combined effect of relic density and DD constraints, this model is ruled out way beyond the region favoured by the flavour anomalies. One could possibly overcome this problem by assuming a non-standard cosmological history of the early universe, see e.g. [106,107,[123][124][125], such that the DM abundance is diluted to the extent that the region compatible with the flavour anomalies in the last plot of Figure 11 becomes partly viable.

S IIB , Majorana singlet DM
Here, we consider a variant of the previous model featuring Majorana rather than Dirac DM. As shown in Figure 12, this time moving from Dirac to Majorana DM does not open new viable regions of the parameter space. The Majorana nature of the DM particle eliminates the Z-penguin contribution to the scattering with nuclei discussed in the previous subsection. This noticeably relaxes the DD constraints for the case of large couplings to the top as in our scenario i ), shown in the first line of Figure 12. However, irreducible QCD contributions to DD still impact the parameter space even for Majorana DM. Furthermore, the region compatible with the flavour constraints typically correspond to overabundance of DM. In the Majorana case, the DM annihilation cross-section is even more suppressed, because of velocity dependence, with respect to the case of Dirac DM. As a result, this model is still not viable unless a non-standard cosmology provides additional DM dilution. As already pointed out, this outcome follows from to the fact that the DM is coupled only with one NP state. This is the main difference with e.g. model F IA; 0 that, as we showed in Section 4.2, easily fulfills all constraints in the case of Majorana DM.

S IIIA; -1/2 , Majorana triplet DM
To conclude our overview of scenarios with distinct phenomenology, we illustrate two models with DM belonging to a SU (2) L triplet. In fact notice that, for all models in Tables 2 and 3 featuring a complete singlet, gauge invariance allows to substitute to the singlet an SU (2) L triplet with zero hypercharge. As we will see below, the consequent change in the phenomenology of our models is dramatic.
We start with the model S IIIA; -1/2 , whose DM candidate is part of the fermion triplet Ψ L , (1, 3, 0). The field Ψ Q also transforms as an SU (2) L triplet, as well as a colour triplet: (3, 3, 2/3). The mediator is a complex scalar doublet (1, 2, −1/2). The Lagrangian reads: As mentioned above, the DM candidate is the neutral component of the Majorana field Ψ L . The combined constraints on this model are shown in Figure 13. We consider again the single assignment M Ψ Q = 1400 GeV, as our usual benchmark values 700, 1100 GeV are mostly excluded by recasting the LHC searches in Refs. [74,75] in terms of the process pp → Ψ Q Ψ Q → qq + ΦΦ → qq + µ + µ − + / E T . The S IIIA; -1/2 model features the weakest correlation among flavour/LHC and DM observables. As we can see from the Lagrangian, the DM is coupled only with the colour singlet Φ. Being in addition a Majorana fermion, the only contribution to SI interactions comes from loop diagrams involving the charged components of the DM multiplet Ψ L as well as the W, Z bosons [58,116]. This kind of interactions lead to cross-sections still below current experimental sensitivity [126,127]. For this reason no DD exclusion region appears in Figure 13. For what concerns the relic density, belonging the DM candidate to a triplet, it features a very efficient and possibly Sommerfeld-enhanced annihilation cross-section into gauge boson pairs (cf. e.g. the third diagram in Figure 2) such that the CMB bound Ω DM h 2 0.12 is saturated only Figure 13: Summary of the constraints for the model S IIIA; -1/2 with triplet Majorana DM. We have considered the values 1, 2, 3 for the |Γ Q µ | coupling. Contrary to the other models, DM constraints do not depend on the individual values of the Γ Q s and Γ Q b couplings, hence have not made the usual distinction between the configurations i ) and ii ). The colour scheme is as defined in Figure 5. In addition, the yellow region is excluded by LHC searches of disappearing tracks, see the text for details.
for DM masses of the order of 3 TeV, far from the region compatible with the fit of flavour observables. In the regime shown in the plots, the DM is always underabundant (unless some non-thermal production mechanism is assumed), irrespective of the values of the masses and couplings of the NP fields and, hence, no relic density exclusion appears. From Figure 13, we can also see that the region excluded by LHC searches for events featuring missing energy (in this case pp → ΦΦ → µ + µ − + / E T [71]) is not very pronounced. However, compared to the previous models, these plots feature a new type of excluded region, filled in yellow. This bound corresponds to the negative results from LHC searches [78,79] for the disappearing charged tracks that, in the model under consideration, would be associated to the pair production of the electrically charged components of the electroweak multiplet the DM belongs to, see the related discussion in Section 3.2. In the minimal setup considered in this work, the different states composing the DM multiplet have loop-suppressed O(100) MeV mass splitting determined by electroweak gauge interactions. As a consequence, the charged DM partner is long-lived and decays into final state particles which are too soft to be detected or would, hence leading to disappearing track events rather than prompt jets or leptons and missing energy events. For the S IIIA; -1/2 model we can directly apply the disappearing-track bound obtained for the case of a supersymmetric Wino [79] which translates into a limit of the triplet mass M Ψ L ≥ 490 GeV. As we can see, this latter bound completely covers the region fitting the flavour anomalies.

F IIIA; -1/2 , Real-scalar triplet DM
The last model we consider is analogous to the previous one with the role of fermion and scalar fields reversed. It is described by the following Lagrangian:  triplet Φ L , and the mediator is the fermion doublet Ψ.
The effect of the combined constraints on the F IIIA; -1/2 model is shown in Figure 14. For analogous reasons to those illustrated in the previous subsection, also in the case of scalar triplet DM, we notice the absence of bounds coming from DM direct detection and relic density. Detailed studies of the DM phenomenology of real scalar triplets have been conducted e.g. in Refs. [128,129]. A notable difference with respect to the previous model S IIIA; -1/2 emerges, on the contrary, for what concerns LHC bounds. Indeed, the bounds from missing energy events (specifically on pp → ΨΨ → µ + µ − + / E T [71]) impact a larger (orange) region of the parameter space, compared to the analogous model with fermionic DM. This is again due to the fact that the Drell-Yan production cross-section of the fermion pair Ψ is substantially larger than the one of a scalar pair with the same mass and the same quantum numbers under the SM gauge group. On the contrary the bound from disappearing tracks (here we show the limit as recasted for the case of a scalar triplet in Ref. [128]) is weaker for scalar DM, again due to the different production cross section. As a consequence, we notice the presence of (narrow) regions of the parameter space compatible with the flavour anomalies, provided that |Γ L µ | 3. We expect that these unconstrained regions can be tested employing future LHC data by a combination of searches for disappearing tracks and searches for events with soft leptons and missing energy like those in Refs [72,73].

Summary and conclusions
In this work, we have presented a systematic study of minimal scenarios providing a viable fit of the observed anomalies in semileptonic B-meson decays and simultaneously solving the DM puzzle thanks to a particle candidate that can achieve, through the thermal freeze-out mechanism, a relic density compatible with the experimental determination from CMB anisotropies. In this minimal setup, the SM spectrum is extended by three new states, either two scalars and a fermion, or two fermions and a scalar, coupled, according to gauge invariance, with left-handed muons and quarks of the second and third generation. All these new fields, including the DM candidate, are present in the loop diagrams associated to the NP contributions to the rates of B-meson decays, as shown in Figure 1. This kind of models hence features an interesting connection between flavour and DM physics. To our knowledge, the present work shows for the first time a complete classification of the possible models of this kind which can be elaborated, depending on the quantum numbers of the new fields. The details of the considered setup are given in Section 2 and the outcome of such a classification is summarised in Tables 2 and 3. In Section 4, we have studied in detail a selection of these models encompassing a large variety of scenarios. Among the models we chose, four possible natures for the DM candidate (namely real scalar DM, complex scalar DM, Dirac and Majorana fermionic DM) are represented. Furthermore, our selection includes examples with the DM field being a singlet of SU (2) L , as well as cases of DM belonging to an SU (2) L doublet or triplet. Following the strategy described in Section 3, for each model we have performed a fit to the B-physics anomalies and used the results to define benchmark assignments for the couplings of the new particles with quarks as well as the mass of one of the (non-DM) NP fields. We have then studied, in terms of the remaining parameters, a broad range of constraints: bounds from searches for the new states at the LHC, DM relic density, DM direct detection and, when appropriate, DM indirect detection.
The results of this analysis have been presented in detail in Section 4, and the general lessons that we can extract from it can be summarised as follows.
• A good fit to the flavour anomalies is possible if the product of the couplings of the NP fields to bottom and strange quarks is moderate Γ Q ∼ 0.15 (larger values would be in conflict with constraints from B s mixing) and consequently the coupling to muons must be rather large Γ L µ 2 (to fit the anomalies at the 2σ level, cf. Section 3.1). This has important consequences for DM phenomenology: if DM belongs to one of the two fields coupled to muons, annihilations into muons are very efficient in depleting the DM abundance to (or below) the observed value. Moreover, electroweak penguin diagrams like the one depicted in Figure 4 can give a large contribution to the DM-nucleon scattering cross section relevant for DM direct detection.
• As a consequence, in the cases with DM coupling to muons, especially if it belongs to the field (Φ or Ψ) that acts as "flavour messengers" in Figure 1, we observe a high degree of correlation among our observables, namely the couplings of the NP fields to SM fermions simultaneously control DM, flavour, and collider observables (this is the case for instance of models F IA; 0 and S IA featuring singlet DM). Furthermore, the relic density constraints can be easily satisfied in the region of the parameter space that fits the flavour anomalies.
• However, strong constraints from DM direct detection would substantially rule out these scenarios, unless DM is a Majorana fermion or a real scalar (cases for which the most relevant DM-nucleon operator vanishes) or it is a complex scalar with a mass splitting > O(100) keV between its two components (making the scattering inelastic). In fact, model F IA; 0 with Dirac DM is completely excluded (cf. Figure 5) while it is among the most favourable scenarios from an experimental perspective if DM is Majorana, see Figures 6,7. Similarly, S IA is a viable option and fits well both DM and the flavour anomalies only if the above-mentioned mass splitting is assumed, cf. Figures 9, 10. Interestingly, the viable regions of the parameter space of these models are already partially constrained by LHC searches for jets/muons and missing energy and by direct detection, hence they have good prospects of being tested by next generation detectors like XENONnT [130] and future runs of the LHC.
• For models where DM is still a singlet but couples only to quarks (such as in the example S IIB ), DM annihilation is typically not efficient enough and the fit of the flavour anomalies points toward regions of the parameter space where DM is overproduced. These cases are then typically excluded by the relic density constraint independent of whether DM is a Dirac or a Majorana field, cf. Figure 11 and 12, but they could be possibly viable within modified cosmological histories of the early universe providing additional DM dilution.
• Our analysis also shows that a combined fit of DM and flavour anomalies favours scenarios where DM is a singlet of the SM gauge group. If DM is instead part of an SU (2) L multiplet the correlation between relic density and flavour observables is lost, as DM annihilation mainly proceeds through gauge interactions, thus independently on the couplings with quarks and leptons. Furthermore, one should rely on a non-thermal DM production mechanism since, in the regions of parameter space where a viable fit of the B-anomalies is achieved, the DM is always underabundant in light of its very efficient annihilations into gauge bosons, see models F IB; -1/3 , S IIIA and F IIIA . The main challenge to this kind of models comes from LHC searches. In particular, in case of DM belonging to an SU (2) L triplet the interesting signature of disappearing charged tracks excludes or drastically restricts the regions of the parameter space compatible with the flavour physics anomalies, see Figures 13 and 14.
As mentioned in the introduction, the present exercise did not aim at proposing "realistic" BSM scenarios, rather at highlighting the minimal ingredients that a more fundamental theory may need to include if the new physics (possibly) behind the B-physics anomalies is indeed related to the DM sector. The above analysis studied the role and the phenomenological consequences of such minimal building blocks. These scenarios could be easily extended to include more particles and interactions. In particular, additional vectorlike fermions or scalars, mixing through a SM Higgs vev with the fields considered here, would also induce operators involving right-currents that may provide an even better fit to the b → s data. Similarly, this would introduce couplings to right-handed muons that can realise chirally-enhanced contributions to the muon g − 2 and thus a natural fit of the observed anomaly, see e.g. [68]. Moreover, a more realistic flavour structure of the couplings (rather than our ad hoc assignment) could be considered, possibly following from some flavour symmetry or other models explaining the observed hierarchies of fermion masses and mixing. Within frameworks of such kind, one could find correlations between our observables and flavour processes in other sectors (e.g. s − d transitions), and thus additional constraints and handles to test the scenarios we considered.