Fermionic Singlet Dark Matter in One-Loop Solutions to the $R_K$ Anomaly: A Systematic Study

We study the dark matter phenomenology of Standard Model extensions addressing the reported anomaly in the $R_K$ observable at one-loop. The article covers the case of fermionic singlet DM coupling leptophilically, quarkphilically or amphiphilically to the SM. The setup utilizes a large coupling of the new particle content to the second lepton generation to explain the $R_K$ anomaly, which in return tends to diminish the dark matter relic density. Further, dark matter direct detection experiments provide stringent bounds even in cases where the dark matter candidate only contributes a small fraction of the observed dark matter energy density. In fact, direct detection rules out all considered models as an explanation for the $R_K$ anomaly in the case of Dirac dark matter. Conversely, for Majorana dark matter, the $R_K$ anomaly can be addressed in agreement with direct detection in coannihilation scenarios. For leptophilic dark matter this region only exists for $M_\text{DM} \lesssim 1000 \, \mathrm{GeV}$ and dark matter is underabundant. Quarkphilic and amphiphilic scenarios even provide narrow regions of parameter space where the observed relic density can be reproduced while offering an explanation to $R_K$ in agreement with direct detection experiments.

Our work builds on a class of models introduced in [13] to address the R K anomaly. The authors analyze a scenario where contributions to b → sµ µ are realized at one-loop level via three new particles charged under the SM gauge groups. Such new physics scenarios can also provide sizable contribution to the muon magnetic moment ((g − 2) µ ). Allowing for representations up to the adjoint, they find and analyze 48 different models. Several of those models include an electrically neutral color singlet, which, if stabilized by a symmetry, can be a valid DM candidate. In this work, we enlarge the analysis of this class of models and address the question if both the R K anomaly and the observed DM relic density can be explained within this framework for models including a fermionic singlet. The connection of this class of models and DM has been addressed previously, for instance in [14][15][16][17][18]. All of these specify on one out of the 48 available models. References [14,16] discuss a scalar singlet DM realization (bIA model), while [15] discusses DM in the context of a leptophilic fermionic singlet setup (bIIA) for Majorana DM. The article [17] discusses the bIIA setup in the context of scalar doublet DM. Reference [18] addresses fermionic Majorana DM in the aIA model specification. The connection of B-anomalies and DM has also been studied in different one-loop setups [19] and tree-level realizations [20,21].
Typically, this setup requires a relatively large coupling of the new particles to the second lepton generation to explain the R K anomaly, since, in comparison to tree-level realizations, the one-loop suppression has to be compensated and, additionally, the coupling to quarks is constrained by B 0 -B 0 mixing. This in turn can lead to an underpopulated dark sector, as the strong coupling to the second lepton generation tends to delay the DM freeze-out. Moreover, this coupling can generate a large vector current coupling of DM to the Z-boson at one-loop level that is tightly constrained from the XENON experiment [5]. However, in for instance models with a fermionic singlet DM candidate, DM can be either a Dirac or Majorana fermion, where the latter has a vanishing vector coupling to the Z-boson and weakens the constraints from B 0 -B 0 mixing. Thus, these models seem especially promising and are the main subject of our studies.
This article is structured as follows: In Section II, we review the model and constraints presented in [13]. Subsequently, we discuss the implications of the models on the relic density and direct detection experiments in Section III, where we also provide a summary of our analysis strategy for our numerical analysis. The results of this analysis for each model are presented in Section IV. In Section V, we conclude.

A. Setup and Classification
The model classes proposed by [13] as a one-loop solution to the b → sµ − µ + anomalies are distinct in their particle content. In either of the two model realizations, three additional particles are added to the SM. In realization a, two heavy scalars, φ Q and φ L , and a vector-like fermion ψ are present, whereas in realization b there exist two vector-like fermions, ψ Q and ψ L , and a heavy scalar φ . The indices Q and L denote their coupling to quarks/leptons respectively. The corresponding Lagrangians for each realization are L a int = Γ Q iQ i P R ψφ Q + Γ L iL i P R ψφ L + h.c. (1) where Q i and L i denote the left-handed quark/lepton doublets, while i is a flavor index. The Eqns. (1) and (2) reveal that the additional heavy particles couple only to left-handed SM fermions. This is a phenomenologically driven feature, implemented to ensure that C 9 = −C 10 holds, which is the preferred scenario to fit the B anomalies, where C 9 and C 10 are Wilson coefficients of O 9 = (sγ ν P L b)(μγ ν µ), O 10 = (sγ ν P L b)(μγ ν γ 5 µ) .
The model classification up to the adjoint representation can be extracted from Table I In this context, X is defined as the hypercharge of ψ in model class a, while it is defined as the negative hypercharge of φ in model realization b. The parameter X can be freely chosen in units of 1 /6 in the interval X ∈ (−1, 1).
In this work, we want to analyze whether the model classes presented above contain a viable DM candidate, while still being able to explain the B anomalies. A dark matter candidate is constrained to be a colorless, electrically neutral, massive particle 1 . This statement alone eliminates one half of the 48 possible model configurations, since in the categories C and D there is not a single colorless particle. The remaining 24 models can be classified in terms of the properties of their DM candidates. In this article, we limit ourselves to models with a fermionic singlet dark matter candidate, which amounts to five models where the 1 In fact, DM could be colored and exist in form of eventually colorless bound states. We, however, only consider single particle DM. Please note also, that in principle DM could possess a small electric charge, such as in scenarios of millicharged DM [22,23]. Table I. All possible choices for the combinations of representation of the new particles such that they allow for an one-loop contribution to b → sµ + µ − . The upper sign of ± belongs to a-type models, the lower to b-type models.
DM can be either Majorana or Dirac fermion. All singlet DM models are categorized in

B. Constraints on new Yukawa Couplings
Following [13], we obtain constraints on the couplings Γ µ from constraints on the Wilson coefficients C 9 and C BB , which are obtained from global fits of LFUV observables and B-B-mixing respectively, which is generated by the effective operator O BB = (s α γ µ P L b α )(s β γ µ P L b β ) .
The Wilson coefficients read C box, a 9 = −C box, a where The (BB) can be extracted from Table III. As in [13], we neglect the influence of photon penguin diagram contributions to C 9 and therefore assume C 9 ≈ C box 9 . The 2σ bounds on C 9 = −C 10 and C BB are [24,25] Starting from these premises, we can construct an upper bound on Γ s Γ * b from B-B mixing and use this to construct a lower bound on Γ µ by taking into account the bounds on C 9 .
We assume mass degeneracy between the new non-dark matter particles φ Q and φ L in a-type models and ψ Q /ψ L and φ in b-type models respectively. 2 For convenience, we further introduce the dimensionless parameter κ, which we define as quantifying the mass gap between DM and the non-DM exotic particles in the model.
Additionally, we restrict our analysis to single-component DM scenarios. This translates into a lower bound on the decay rate of the heavier dark sector particles to ensure their decay proceeds sufficiently fast. For the mass splitting, this in turn translates to (κ − 1)m DM > m π [26]. Within this model it is possible to construct a solution to the anomalous magnetic moment of the muon, commonly dubbed (g − 2) µ . In these models we have a contribution to a µ = (g−2) µ/2 as where the group factorsη a µ /η a µ and χ a µ can be extracted from Table III and the functions F 7 (x) andF 7 are characterized as Recently, the g-2 collaboration updated the earlier results of [27] and reported on a 4.2 σ deviation from the SM, where the difference ∆a µ ≡ ∆a Exp µ − ∆a SM µ amounts to [28] ∆a µ = (251 ± 59) · 10 −11 .
This translates to upper and lower bound on the leptonic Yukawa coupling Γ µ summarized in Table VII.
The constraints on the new Yukawa couplings Γ µ ,Γ b and Γ s are summarized in Table VI. In general, Dirac and Majorana versions of the models can lead to different constraints on the Yukawa couplings, 2 We comment on the effects of lifting this assumption for the DM phenomenology at the beginning of section IV since additional contributions to the B-B-mixing and b → s l + l − transitions can arise. These contributions, however, only occur in the case of model aIA 3 because of the constellations of SU (2) representations present in this model. This feature is thus ultimately incorporated in the Wilson coefficients C BB and C box 9 , together with Table III, where the product χ M (BB) η M (BB) is non-zero only in the case aIA. In general, the the constraints of the new Yukawa couplings read The coefficient functions B model i (κ) are model dependent and summarized in Table VI.

C. Constraints on Parameters from the Scalar Potential
The extended scalar sector of the models leads to additional terms in the scalar potential. These parameters are mass parameters µ i as well as quartic couplings λ . In this work, we assume that no other scalar than the SM SU(2)-doublet scalar field H acquires a vacuum expectation value since a non-zero vacuum expectation value in the dark sector would break the stabilizing Z 2 symmetry. The resulting condition on the model parameters is given specifically below.
Further, we demand perturbativity from all allowed quartic couplings. This yields We also want to constrain ourselves to scenarios, where the vacuum is stable at tree-level. Since the particle content in the scalar sector differs between model realizations a or b, the implications for the scalar quartic couplings differ as well.

a-type Models
There exist two other scalar doublets φ Q and φ L in a-type models in addition to the SM SU(2)-doublet H. The most general scalar potential, which is invariant under the DM stabilizing symmetry in this scenario is consequently Adopting the limits from [29], the vacuum stability bounds at tree-level for this kind of potential read where α ∈ [φ Q , φ L , H] and i ∈ [1, 2, 3]. Since we assume mass degeneracy within the non-DM dark sector, the Higgs portal couplings λ φ L ,H,2 and λ φ Q ,H,2 vanish. Another important feature of the model is the DM stabilizing symmetry Z 2 , under which all dark sector particles are oddly charged. For this symmetry to be intact, we require the other Higgs portal couplings to satisfy

b-type Models
In addition to H, there exists another scalar doublet φ in b-type models. The most general form (again respecting the DM stabilizing symmetry) of the scalar potential in this scenario is therefore For this kind of models we adopt the limits of [30,31], leading to As in a-type models, we can infer constraints for the Higgs portal couplings from mass degeneracy of the non-DM exotic particles and the requirement that the new scalar does not acquire a non-zero vev, yielding

III. DARK MATTER PHENOMENOLOGY
In this section, we analyze whether the aforementioned models can simultaneously create the correct DM relic density Ω DM h 2 = 0.120 ± 0.001 (at 68% CL) observed by Planck [1] and withstand the bounds set by DM direct detection (DD) experiments like XENON [4]. To carry out our phenomenological analysis, we use FEYNRULES [32,33] and its interface with FEYNARTS [34] and FORMCALC [35] and/or FEYNCALC to calculate contributions to direct detection cross sections up to one-loop level. We also use another FEYNRULES interface with CALCHEP [36] to compute direct detection cross sections with MI-CROMEGAS, which uses said CALCHEP output. We additionally use MICROMEGAS to numerically solve Boltzmann's equations [37,38] and obtain the relic density. Since the Yukawa coupling required to generate the observed R K is sizable, the DM production mechanism present in this model can only be the standard freeze-out mechanism.

A. Relic Density
The relic density can be estimated by the effective annihilation cross section σ eff of DM. Depending on the parameter κ, which parametrizes the mass ratio of the various particles in the dark sector according to  [39,40] Here, we use x = m DM /T and M Pl = 1.2 · 10 19 GeV is the Planck mass. The yieldỸ is defined as the sum over all coannihilating particlesỸ while the Yield Y i of the particle species i is related to its number density n i via Y i = n i /s, with the entropy density s of The quantity g * ,eff combines the entropy degrees of freedom g * s and the energy degrees of freedom g * via Since we assume mass-degenerate unstable dark sector particles with a mass of κM DM the yields can be expressed as Here, g DM and g non-DM refer to the degrees of freedom of the DM and the non-DM dark sector particles, respectively. For instance, in a b-type model with ψ L DM we have g DM = g ψ L and g non-DM = g ψ Q + g φ . The effective annihilation cross section is given by a weighted sum of the thermally averaged cross sections of the various annihilation channels where i and j run over all dark sector particles and σ i j v i j is the thermally averaged cross section for for the annihilation of the dark sector particles i and j. 4 Note that MICROMEGAS 5.0 does not validate that conversions between the dark sector particles are efficient but assumes them to be efficient enough to allow for the treatment described above. We numerically verified the efficiency of 2 ↔ 2 conversion processes and 1 ↔ 2 (inverse) decays between dark sector particles for coannihilating scenarios κ 1.2 during freeze-out, i.e.
Γ conversion > H for T M DM /30. Reference [41] discusses the efficiency of conversions in a slightly simpler but similar setup and find that effects of out-of equilibrium conversions only arise for couplings Γ i 10 −6 , which is far below the couplings considered in this work. Due to the singlet nature of DM, the direct DM annihilation proceeds via the Yukawa couplings Γ i connecting DM to the SM leptons and/or quarks 5 , while the heavier dark sector particles may be charged under the SM gauge groups and therefore can additionally annihilate via gauge interactions. Various annihilation channels are depicted in Figure 1. The relic density of a particle species 6 can be estimated from the leading order velocity contribution of the annihilation cross section σ eff v ∼ σ eff,0 T n [42] where x f ∼ 25 parametrizes the freeze-out temperature, m B is the typical baryon mass and Ω B and Y B describe the energy density and the yield of baryons today, respectively. Assuming that the annihilation cross section is dominated by the annihilations of dark matter into SM leptons and quarks of negligible mass 7 via the new Yukawa couplings, the leading order contribution of the thermally averaged cross section results in Further, we find n = 0(n = 1), corresponding to s-wave(p-wave) annihilation, for Dirac (quarks) exclusively. 6 In case of Dirac DM the total relic density is given as the sum of the particle and antiparticle contribution. 7 This scenario arises frequently in models with ψ L , ψ-DM or for ψ Q -DM significantly heavier than the top. coefficients are equal to 1.
Given the σ −1 eff,0 scaling of the relic density in Eq. 36, we expect to generate the observed relic density in such a scenario neglecting the logarithmic mass dependence of x f for In the following we derive some implications for limiting behaviors of this condition for the various models.

Leptophilic DM (bIIA,bVA)
In this scenario Eq. (38) results in a scaling behavior Comparing this result with the lower bound on Γ µ for a successful explanation of the R K anomaly, given in Eq. (19), we find a lower bound on κ > κ m 0 for a possible simultaneous explanation of R K and DM of

Quarkphilic DM (bIIB,bVIB)
In this scenario the observed relic density is reproduced for The product of Γ s and Γ b is constrained from above by B −B-mixing, given in Eq. (18). The annihilation cross section is given by the sum Γ 2 s + Γ 2 b . Thus, it cannot be constrained from B-B mixing since either the second or third generation coupling could be arbitrarily small. However, we obtain an upper bound on the annihilation cross section by setting Γ s ∼ 0 and Γ b = 4π at its perturbative limit 8 For the four different configurations we find

Amphiphilic DM (aIA)
Since both, couplings to leptons an quarks, are present in this scenario, we can apply the limits derived in the two scenarios above if Γ µ Γ q or Γ q Γ µ .

B. Direct Detection
In this section we discuss the direct detection limits arising for the exemplary case of a fermionic singlet DM particle χ interacting with the SM quarks via a vector (V µ ) or scalar (A) mediator. Dark matter direct 8 Note that the second generation coupling Γ s cannot be close to its perturbative limit, as it is constrained from D-D mixing. instance XENON1T, LUX, PICO and also IceCube put bounds on the DM-nucleon cross section, we shortly discuss the estimates of the spin-independent (SI) as well as spin-dependent (SD) parts of this quantity for the different diagram structures present in this work. For this discussion, we mainly use the results of [43], [44] and [45]. We put the main focus on dependence of the cross sections on the DM mass and the couplings.
An important note is that the bounds presented by the aforementioned experiments are given under the assumption that the total DM relic density is constituted of the particle in question. In the case where a model underproduces the relic density of a DM candidate, the DD bounds have to be rescaled such that only the produced fraction of the observed relic density Ω DM h 2 is regarded. Thus, the relation must hold for a model to be in accordance with DD bounds. Bounds for underproduced DM are relaxed, as the allowed DM-nucleon cross section rises by a factor Ω DM/Ω χ > 1. If the RD is in turn overprocuded, the bounds are tightened in our setup. Note however, that this part of the parameter space is already excluded by the RD and therefore DD results in these regions are of no special interest.
We estimate the results for t-channel DD interactions 9 as shown in Figure 2, where an EFT reduction is made, which is eventually matched to the operator structure in Eq. (42). If the DM-quark operator structure differs from this arrangement, as for instance in an s-channel diagram, this operator structure must be translated to the t-channel structure via Fierz transformations. This is done because for the t-channel operators, DM and SM particle content strictly factorizes.
In this work, we evaluate the interaction within the blobs in Fig. 2 up to one-loop level with FEYNARTS.
The effectiveχ χH/Z/γ/g -vertices are then implemented in the CALCHEP-files, which are used by MI-CROMEGAS to calculate the DM-nucleon cross section.
Typical t-channel operators with bosonic mediators contributing to DD are structured as where the Π med is the propagator of the mediator and O N denotes the nucleon-level operator. As we are interested in the nucleon rather than the partons, a summation over the quark content is necessary. This procedure, even up to the level of the whole nucleus, is explained in detail in the appendix of [43] and references therein.

Fermionic DM χ and Scalar Mediator A
In the fermionic DM and scalar mediator case, the Lagrangian yields where the factor of 1 2 enters in the Majorana case. In this type of interaction, the combinations scalar-scalar (s, s),scalar-pseudoscalar (s, p), pseudoscalar-scalar (p, s) and pseudoscalar-pseudoscalar (p, p) between the DM and the nucleon operators O DM and O N arise. The DM-nucleon cross sections for each combination scale like where µ χN = m χ m N /m χ +m N denotes the reduced mass of DM χ and the nucleus N . Interactions that are scalar on the nucleon operator contribute to the SI cross section, whereas SD contributions arise from interactions that are pseudoscalar on the nucleon operator side. Note here especially that the (s, s) interaction is the only interaction unsuppressed by the relative velocity v. The mixed terms (s, p) and (p, s) are suppressed by two powers of v, while the (p, p) interaction features a suppression of v 4 . In this work, those interactions are mainly SM-Higgs-exchange induced.

Fermionic DM χ and Vector Mediator V µ
In the fermionic DM and vector-boson mediator case, the Lagrangian reads where the factor 1 2 enters only in the Majorana case. A special caveat here is that g χv vanishes in the Majorana case, since terms of the formχγ µ χ are forbidden. For the possible interactions, vector-vector (v, v), vector-axial vector (v, a), axial vector-vector (a, v) and axial vector-axial vector (a, a), we obtain Analogously to the scalar mediator case, the nucleon current containing γ 5 , the axial vector current, leads to SD contributions, while the SI contribution stems from the vector current. The velocity suppression in these configurations is, however, not quite analogous. The (v, v) and (a, a) interactions are unsuppressed, while the mixed terms (a, v) and (v, a) are suppressed by two powers of v. Therefore, there does exists a contribution to the SD cross section unsuppressed by velocity.

s-channel Lagrangians and Fierz Identities
In the case of s-channel-type direct detection, it is possible to relate the Lagrangian to the t-channel interactions whose cross sections we presented in Section III B 1 and III B 2 via Fierz identities. The schannel type DD interactions that are typically present in quarkphilic models are mediated by a scalar.
Starting from the structure of the Lagrangian we obtain the effective Lagrangian [44] by integrating out the mediator. Further manipulation with Fierz transformations yields where velocity suppressed combinations are neglected. Eq. (50) suggests that in the limit |λ s | = |λ p |, the scalar-scalar and tensor-tensor operators vanish completely. The remaining quadrilinears can be related to the scaling of the cross sections presented in Sections III B 2 and III B 1. Note that both vector-vector and tensor-tensor contributions vanish in the Majorana case.

Twist-2 Contributions to the SI DM-Nucleon Cross Section
In some Majorana versions of models presented in this work, there is no unsuppressed contribution to the SI cross section from quark operators at leading order. In this case, twist-2 operators become important.
While the twist-2 quark operators come from a higher order propagator expansion of the tree-level quark diagrams, also gluonic twist-2 operators generated by the box-diagrams shown in Fig. 3 have to be considered. These diagrams are taken into account by MICROMEGAS automatically. In this case, the SI cross sections receive additional contributions from these diagrams, which scale like [45] We will refer to these analytic expressions in Section IV, where we present the numerical results, as they play an important role in the SI direct detection procedure in quarkphilic DM Majorana models. Note that in the expressions above, m med corresponds to the mass of the colored particle directly coupling to DM.

C. Indirect Detection
In this section we shortly discuss indirect detection limits for the model classes presented in this work.
Majorana type models are typically unconstrained by standard indirect searches due to the kinematic suppression of the pair annihilation cross section σ (χ χ →f f ) ∼ v 2 . However, this suppression is lifted for dark matter annihilations into a three-body final state containing a fermion pair and a photonχ χ →f f γ [46]. Such processes, commonly called Virtual Internal Bremsstrahlung (VIB) lead to a special spectral feature that is distinct from astrophysical backgrounds. Analyses searching for these features by FermiLAT [47] and HESS [48] can be used constrain models that allow a DM coupling to light fermions [19,49], such as the model classes discussed in this article.
We calculate the VIB contribution to the annihilation cross section using MICROMEGAS 5.0 and compare the results to the the 95% C.L. upper limits on the annihilation cross section σ v f f γ + 2 σ v γγ for κ = 1.1, 1.01 obtained by [49]. As in the case of DD, we rescale the experimental constraints with a fraction of the observed relic density so that Note that for κ < 2, the contribution of the two-photon final state is negligible. Moreover, for large κ VIB is suppressed by κ −4 and thus we do not present limits for these cases.
Our analysis shows that even for coannihilating scenarios the parameter space of the models analyzed in this work cannot be constrained by VIB. Principally, for large couplings to leptons (in the case of leptophilic/amphiphilic DM) or quarks (quarkphilic/amphiphilic) the DM annihilation cross section can be relatively large. Nevertheless, the corresponding points in the parameter space typically feature a small relic density, which results in relaxed ID constraints because of the quadratic rescaling given in Eq. (52).
Since our numerical results indicate that the relevant parameter space for Dirac models is excluded by DD already, we do not apply the standard ID searches constraining the DM annihilation into a pair of SM fermions for these models.

IV. RESULTS OF THE NUMERICAL ANALYSIS
In this section we present the results of our numerical analysis. The section itself is further divided into separate discussion of the individual models. Since the discussed models all contain a singlet fermion DM • We restrict ourselves to parameter sets respecting the DM stabilizing symmetry and satisfying the vacuum stability conditions as well as perturbativity bounds (see Sections II B and II C). The bounds on R K and (g − 2) µ , however, are not regarded as stringent bounds on our parameter sets as we perform a scan over Γ µ .
• We keep the product Γ s Γ * b at the upper bound 10 . This generally facilitates the realization of a solution to the R K anomaly, as smaller leptonic couplings Γ µ suffice in this case. As only the product is fixed, we differentiate between two choices of coupling structures for the new quark Yukawa couplings: democratic and hierarchical. In the democratic scenario we choose the second and third generation couplings Γ s and Γ b to be almost identical 11 In the hierarchical setup we impose a strong flavor structure. As Γ b is unconstrained, it is possible to set its value up to the perturbative limit. In the upcoming plots, we present the most extreme hierarchical case of Γ b = 4π ∧ Γ s = B model bs (κ) · M DM/4π GeV to illustrate the possible impact of a nondemocratic flavor structure. 10 If the bound on Γ µ posed by R K is surpassed, Γ s Γ * b is scaled down accordingly to still be able to solve the R K anomaly. 11 By almost identical we mean as close by as the bound on Γ s from D-D mixing allows. Typically, the difference is within O(1%).
• We fix all non-vanishing Higgs portal couplings to 0.1, since they do not alter the DM phenomenology in a substantial manner.
• We assume mass degeneracy between the new non-dark matter particles φ Q and φ L in a-type models and ψ Q /ψ L and φ in b-type models respectively.
In addition to the scenarios presented in this article, we have studied an examplatory case for each model class (lepto/quark/amphiphilic) allowing for a non-degenerate, non-DM dark sector. For quarkphilic and leptophilic scenarios, only coannihilating scenarios are affected by a non-degenerate, non-DM dark sector, while amphiphilic models can be affected in both setups. Adopting the point of view that the mass degenerate setup refers to the case of the smallest possible dark sector masses involved for a given κ, all regions of viable parameter space are shifted to smaller DM masses when considering a non-mass-degenerate, non-DM dark sector. While an analysis of the cases specified above could be interesting, a detailed study of this scenario is beyond the scope of this paper.
In all plots, the area allowed by direct detection (DD), R K , and muon g − 2 (∆a µ ) are presented in colors blue, gray, and green respectively. The lines belonging to DM relic density (RD) correspond to the observed relic density Ω DM h 2 = 0.120 ± 0.001 in the universe. On the one hand, the parameter space to the left of these lines features underproduction of relic density, while on the other hand the region to the right of the line features overproduction of DM.

A. bIIA
Following Table I, the representations of the dark sector particles in bIIA are As ψ L qualifies as the DM candidate, bIIA belongs to the class of leptophilic DM models.

Dirac DM
The numerical results obtained for bIIA Dirac are summarized in Figure 4. Figure 4 indicates that of the depicted mass configurations only κ = 15 principally allows for a solution of R K , since the corresponding  Coannihilating scenarios are not always dominated by the aforementioned annihilations and therefore do not behave as plainly as their non-coannihilating counterparts in this model. As can be seen in Fig.   4, there is a lower mass threshold on successful RD reconstruction stemming from annihilations of ψ Q via the strong gauge coupling or via the new quark Yukawa and also φ via the Higgs portal (see Figure   1). The threshold exists also in the case of vanishing Higgs portal coupling and the new quark Yukawa couplings Γ s and Γ b , since the contribution from the strong gauge coupling g 3 always exists. Typically, Higgs portal interactions are dominated by the other two types of interaction and therefore we choose to fix the corresponding coupling λ φ ,H,1 during the entire scan. The visible κ-dependence of the threshold is explained by the typical coannihilation suppression factor of ∼ exp (−2x f κ) (see Eqs. (33)-(35)): As κ increases, coannihilation channels become more and more suppressed and thus the thermally averaged total cross section shrinks, leading to an increased DM relic density. This way, the threshold shifts towards lower masses for increasing κ.
Furthermore, the annihilation cross section into quarks scales like ∼ Γ 4 b + Γ 2 b Γ 2 s + Γ 4 s at any given κ, which explains the shift of the threshold to higher masses in the hierarchical scenario, since this factor is maximized for large values of Γ b realized in the hierarchical coupling structure. As the masses grow larger, coannihilating scenarios with the correct relic density become dominated by direct annihilations into leptons, as a larger Γ µ is needed to achieve the correct thermally averaged cross section, leading to the same scaling behavior as in non-coannihilating scenarios.
In the democratic setup, another effect can be observed in the coannihilating scenarios: The intersection of the RD lines of κ = 1.1 and κ = 1.01. This effect cannot be understood from annihilations described in Eq. (37) alone, but by taking into account effective conversions of the dark sector particles. We provide a general estimate of the thermally averaged cross section in Appendix A. In there, we also shortly the discuss the κ interval, where such intersections occur.
Direct detection is mediated at leading order at one-loop level in this model, as there are no tree-level vertices with quarks for leptophilic DM. The corresponding leading order diagrams are depicted in Figure   5. It is important to note at this point that the box diagram exhibits one additional heavy dark sector particle propagator. Different from the penguin diagrams also, the final state quarks are of second or third generation, as Γ d = 0, which suppresses the diagram even further due to the reduced fractions of higher generation quark parton distribution functions (PDFs).
Taking DD results into account, we can state that a solution of the R K anomaly is impossible in this setup, since the allowed regions for DD and R K are completely disjoint in both hierarchical and democratic implementation. It is important to stress that DD even excludes configurations where DM is underproduced.
For this model to not be excluded, coannihilations are essential, since only those configurations feature nonoverproduced RD within the allowed DD regions. The tightness of the constraints present in this model is due to the strong effective vector coupling of ψ L to the Z-boson, which in turn lead to a strong constraint from the SI DM-nucleon cross section. This constraint is weaker in coannihilating scenarios by 1 − 2 orders Figure 5. Leading order DM-quark diagrams contributing to the DD cross section of ψ L DM.
of magnitude, which stresses the tension between a solution to R K and DD even more.
In general, two competing effects are at work in this model: On the one hand, the effectiveψ L ψ L Zcoupling becomes smaller with increasing κ, since the masses of the dark sector particles within the loop increase accordingly, thus alleviating the bounds. An increasing κ on the other hand generally increases

Majorana DM
The Majorana version's results are presented in Figure 6. In comparison to the Dirac version, the most striking differences of the Majorana version are the DD bounds, which are weakened by about one order of magnitude in Γ µ . This is mainly due to the vanishing of vector current contribution to theψ L ψ L Z vertex, which is forbidden in the Majorana case. In the absence of this contribution, the most stringent bounds do not come from the SI DD cross section, but from the SD DM-nucleon cross section generated by the axial vector current contribution toψ L ψ L Z vertex, although the SD limits are ∼ 6 orders of magnitude weaker.
The bounds on SD DM-neutron cross section also come from XENON [4], while IceCube puts the more stringent bounds on the DM-proton cross section [50] in this mass range. The DM-neutron cross section tends to constrain the lower mass region up to ∼ 300 GeV, while higher masses are typically constrained by the DM-proton cross section bounds from IceCube. We always apply the strongest bound at any given mass 12 .
As the blue regions allowed by DD now overlap with the gray (R K ) and green ((g − 2) µ ) regions, valid solutions to R K or (g − 2) µ (though not simultaneous) exist. We state here, however, that these solutions do not solve the DM problem, as DM is strongly underproduced in this part of the parameter space. Such constellations thus require a multi-component DM solution 13 that goes beyond this model.
In the democratic setup, the overlap with the R K region exists in the mass region 180 GeV, while a slightly bigger window arises in the hierarchical setup with DM masses ranging up to ∼ 420 GeV for κ = 1.1 12 In the discussion of the upcoming models, the mass regions where the aforementioned experiments are most constraining vary but are not explicitly stated 13 Note that the other DM component(s) must stem from a separate dark sector to the one presented in this work in this case. and ∼ 1000 GeV for κ = 1.01. For a valid (g − 2) µ solution a flavor hierarchy is beneficial, as the allowed mass range for a democratic case is up to less than 100 GeV, whereas in the hierarchical case we observe viable masses up to ∼ 170 GeV for κ = 1.1 and even ∼ 290 GeV for κ = 1.01. The mass regions mentioned above appear to be accessible at the LHC, since they feature a colored particle of a mass less than a TeV.
We briefly review results of searches of similar setups in Section IV F, which suggest such low-massive regions are excluded. A dedicated collider study of this model beyond the contents of Section IV F could be interesting.
Another feature of the Majorana variant is that R K and RD can be reconciled with lower κ values, as Majorana DM annihilations are p-wave at leading order and thus suffer from velocity suppression. This in turn requires a larger Γ µ to achieve to the observed relic density, ultimately pushing relic density lines of lower κ into the gray area. Mathematically quantified, this effect leads to a new lower bound on κ for a reconciliation of R K and RD of κ 4.7 as indicated in Eq. (39).  Table I, the representations of the dark sector particles in bVA are This model can be called the 'sibling' of bIIA, since the fields are almost in the same representations. The only difference is that ψ Q is an SU(2) L -triplet rather than a singlet.  Compared to bIIA, bVA has a greater number of coannihilation channels and therefore features an increased σ v in coannihilation scenarios, leading to a higher mass threshold. As Eq. (39) suggests, the model dependent coefficient B bVA µ (κ) is larger than the one belonging to bIIA. This stems from the fact that there are now also more contributions to B-B mixing and b → sl + l − transitions. This results in a comparably higher κ bVA 0 . Concerning the direct detection results, the main statement from bIIA also applies to bVA. There are no viable solutions to either the R K or (g − 2) µ anomalies in this model. The DD results are virtually the same as in bIIA. This is due to two competing effects. The RD is slightly more depleted for coannihilating scenarios, which weakens the bounds by a factor of Ω /Ω DM , thus rendering the DD bounds weaker than in the bIIA model by O(10%). Figure 8 presents the results for the bVA Majorana version. As indicated in Section IV B 1, there is a strong resemblance between bIIA and bVA. The most interesting difference is the parameter space viable for an R K solution, which shrinks substantially due to the more severe R K bounds. The range is diminished

C. bIIB
Following Table I, the representations of the dark sector particles in bIIB are

Dirac DM
The results of the numerical scans for bIIB Dirac are summarized in Figure 9. As ψ L is the coannihilation partner in quarkphilic DM models like bIIB, Γ µ , the Yukawa coupling directly contributing to R K and ∆a µ is unimportant in scenarios where κ 1.2. As can be seen in Fig. 9, in both democratic and hierarchical versions, κ = 5 and κ = 15 cannot reproduce the observed DM relic density for DM masses > 100 GeV.
This is due to the fact that Γ s/b are not strong enough to deplete the RD by themselves. Further, a Γ µ within the perturbative bounds is not sufficient to overcome the mass suppression by increasing σ v in a way that DM is not overproduced. However, successful relic density reproduction is possible within coannihilation scenarios, where this mass suppression is lowered and thus contributions from ψ L annihilation are sizable enough to offer some viable parameter space. Furthermore, a hierarchical flavor structure between second and third generation quark couplings lowers the influence of Γ µ on the RD, because the thermal freezeout is dominated by direct annihilation into third generation quarks. This behavior is well illustrated by a comparison between Fig. 9 (a) and (b), as no change in the Γ µ -direction is apparent. We also observe the shift of the RD mass threshold to higher masses already described in Sec. IV A 1.
Presumably the most obvious feature of all in this model is the complete absence of parameter space allowed by direct detection. The reasons for this include the aforementioned tree-level diagrams (s and t channel) inducing a sizable contribution to σ SI DD and a vector current contribution to theψ Q ψ Q Z-vertex. All contributions to DD up to one-loop are presented in Figure 10. Note here, that contribution of the effectiveψ Q ψ Q H vertex is small compared to the ones from theψ Q ψ Q Z-vertex. This is the case, although the top-quark Yukawa y t enters in these processes and thus the Higgs exchange is not suppressed by a small Yukawa coupling, which is the case in leptophilic DM models. Moreover, the Higgs portal diagram suffers suppression due to an additional heavy scalar propagator and is thus less important. Figure 11 presents the numerical results of the Majorana DM version of bIIB. The behavior of the relic density is basically the same for Majorana as for Dirac with a slight difference in the mass threshold, which is due to the annihilation being p-wave rather than s-wave.

Majorana DM
With regards to direct detection, the Majorana version exhibits parameter space in agreement with DD limits. As mentioned in earlier sections, the vector current contribution to the effectiveψ Q ψ Q Z vertex vanishes in Majorana DM models, so that the main contribution to the SI DM-nucleon cross section on one-loop level is absent.
The tree-level diagrams (see Fig. 10) are of the scalar mediated s-channel type 14 and can be related to t-channel type DD via Fierz transformations. According to the discussion in Section III B 3, this leads to ) and (a, a) contributions. Since the (v, v) structure does not exist in Majorana models, there is no unsuppressed contribution to the SI DM-nucleon cross section in this case and the twist-2 operators discussed in Section III B 4 become important, which depend on two powers of Γ s/b (see Fig. 3). The mass scaling of σ SI can be derived from Eq. (51). As σ SI effectively diminishes with increasing DM mass and the bound posed by XENON1T softens as ∼ m DM , twist-2 contributions generally constrain lower masses stronger than higher masses. An analogous behavior can be observed for κ. The cross section σ SI is lowered with increasing κ, pushing the threshold of DD to lower masses effectively.
In the democratic scenario, where Γ s/b take moderate values, a considerable amount of parameter space  is left open in coannihilation scenarios. We observe a complementary behavior of SD and SI DD, as SI limits rule out small masses, whereas SD limits rule out larger masses. In the area of small M ψ Q and large Γ µ , RD is significantly depleted by coannihilation channels involving ψ L and thus the bound is softened, while the contributions to the SI cross section do not depend on Γ µ so that high couplings to muons become viable ultimately.
The shape of the SD exclusion can be explained by the DD bound rescaling due to the relic density.
Areas in the parameter space where RD is strongly overproduced feature tighter bounds compared to the underproduced part of the parameter space 15 .
Non-coannihilating scenarios are ruled out by SD DD as the complete parameter space overproduces DM, although they are completely unconstrained by SI DD.
In this model, there is viable parameter space for a solution to R K and DM in the κ = 1.1 scenario for DM masses M ψ Q 1.5 TeV and couplings Γ µ 3. In this window, the correct RD and R K can be For the same reason bVA is considered the 'sibling' model to bIIA, the singlet to triplet shift, bVIB can be considered the 'sibling' of bIIB. 15 Note that this parameter space is ruled out by RD anyways.

Dirac DM
The results for Dirac DM in bVIB are summarized in Figure 12. A comparison between Figs. 9 and 12 shows that the two sibling models bIIB and bVIB expectedly lead to similar results. In bVIB we can observe a more stringent bound on Γ µ for an R K solution than in bIIB. This feature can also be observed in the bIIA vs. bVA comparison, where the triplet-model bVA exhibits more stringent R K bounds. The reason for this is again the increased number of diagrams contributing to the b → sl + l − transitions in a triplet model. The relic densities of the non-coannihilating scenarios κ = 5, 15 are completely unaltered in comparison to bIIB, since only the ψ L gauge representation differ.
Regarding direct detection, it is left to state that the whole parameter space studied in this work is excluded for this model, as is the case in bIIB. The two competing effects, namely rescaling of the DD bounds from lower RD for coannihilation scenarios and more contributions to DD cross section from additional diagrams involving triplet particles can have an O(10%)-effect on the DD bounds (see Section IV B 1). This, however, is not enough to render any of the parameter space viable in the quarkphilic Dirac DM model.

Majorana DM
The results of the Majorana version of bVIB are summarized in Figure 13. This model's results mostly resemble the results obtained for bIIB Majorana, due to the minor differences between the models. The hierarchical scenario is still completely ruled out by direct detection and does therefore not show any significant difference to the hierarchical bIIB Majorana model. Another minor difference is that parameter space allowed by DD in the κ = 5-scenario reaches mass values > 100 GeV, as the mass thresholds of SI DD are generally shifted to higher masses. This is due to stronger underproduction of DM and subsequent relaxation of the DD bounds compared to bIIB Majorana.
Furthermore, the conclusions regarding a solution of (g − 2) µ remain almost the same as in bIIB. There is a possibility for a correct ∆a µ with the downside of underproduced ψ Q -DM. The main difference between the solutions of bIIB and bVIB is the mass range where such a realization is possible. This model shows also that in principle there is a possibility of a solution to (g − 2) µ involving a relatively light ψ Q , as demonstrated in the case of κ = 5. In this case, however, we expect collider searches to become more and more restrictive, as we outline in Section IV F.
As aIA is the only a-type model with a fermionic singlet DM candidate, it makes up it's own category of amphiphilic DM, meaning that DM couples to both quarks and leptons. Thus, aIA possesses with properties of both quark-and leptophilic models. Moreover, this model is the only model that features different R K bounds in Dirac and Majorana versions, as the SU(2) L structure allows for additional diagrams to contribute to the B-B mixing and b → sl + l − transitions (see [13] for a more detailed discussion).

Dirac DM
We present the results for aIA in Figure 14.   Generally speaking, the direct detection results are the most interesting in this model. Since aIA features traits of both quarkphilic and leptophilic models, areas allowed by DD stem from a dynamic interplay between these model characteristics.
As is characteristic for leptophilic models, the DM-nucleon cross section is constrained for large Γ µ , because of the axial vector part of the one-loop diagrams involving leptons.
The shape of the exclusion line in this model, however, differs from leptophilic models. There are two major reasons for this: One the hand, the relic density rescaling is not only influenced by direct annihilation into leptons but also by direct annihilation into quarks in this model. This is the reason for the much more parallel alignment of DD exclusion curves compared to leptophilic models. On the other hand, SD loop-contributions differ because of additional quark loops, which come with an opposite sign compared to leptonic contributions because of their hypercharge structure 17 .
The exclusion from below is a typical characteristic of one-loop contributions involving quarks in the loop, as they are dependent on Γ s/b rather than Γ µ but still sensitive to the masses of the particles running in the loop (and therefore also κ).
The allowed regions visible in Fig. 16 therefore arise due to the interplay of the leptonic and hadronic loop contributions and relic density rescaling of the DD bound.
The above mentioned contributions are all SD but as discussed in the Sections IV C and IV D, SI DD induced by DM-gluon scattering can also occur in quarkphilic models. This effect is visible in coannihilation 16 The R K bound as a function of κ has a global minimum at κ ≈ 1.78 and is divergent at κ = 1, ∞. 17 Note that the additional quark contributions do not necessarily enhance the SD DM-nucleon cross section because of potential cancellations with the leptonic contributions. scenarios, where the allowed regions feature a cut-off at a certain DM mass, such that smaller DM masses are excluded. In the hierarchical scenario, the allowed area for κ = 1.01 disappears completely.
In coannihilation scenarios, we find a 'sweet spot' for a valid DM production in the region M ψ ∈

F. Collider constraints
In this section, we review existing searches of setups similar to the scenarios discussed in the previous sections. While the results cannot be explicitly applied to the setups studied in this work, they can provide an indication of the excluded parameter regions. Note that we do not perform a detailed collider study in this article.
For the leptophilic DM scenarios, i.e. for models bIIA and bVA, ψ Q can be pair-produced at treelevel, subsequently decaying through ψ Q → Qφ † (φ † → ψ LL ). Such channels of ψ L production can be constrained by dilepton+jets+/ E T searches at the LHC [51, 52]: for example, dilepton searches can rule out M ψ L 600 GeV for M ψ Q 800 GeV [15].
For the quarkphilic DM scenarios, such as bIIB and bVIB, ψ Q can be produced in t-channel interactions mediated by the colored scalar φ [53,54]. ATLAS constrained such scenarios through monojet+/ E T searches with an integrated luminosity of 36 f b −1 at √ s = 13 TeV [54]. In our context, for Γ Q = 1, this rules out M ψ Q 600 GeV for M φ ∼ 700 GeV. For lower values of M ψ Q it can rule out even higher values of M φ up to ∼ 1.6 TeV (see Fig. 8 of Ref. [54]). In the model aIA, the DM candidate ψ has tree-level couplings to the SM quarks, leading to the same production channels at the LHC as the bIIB and bVIB models. Thus, one can expect similar constraints on model aIA as well.
In the b-type models, for small enough values of Γ Q or Γ L , colored fermions can become sufficiently long-lived at collider scales and may decay outside the detector. Displaced vertex+/ E T searches at the LHC can be recast into constraints on the mass and lifetime of such a long-lived particle [55]. This can rule out the new physics Yukawa couplings in the range Γ ∼ 10 −2 − 10 −5 for fermion masses 1.8 TeV.

V. SUMMARY AND CONCLUSIONS
In this article, we investigated the DM phenomenology of SM extensions addressing the R K and (g − 2) µ anomalies at one-loop level, which also contain a fermionic singlet DM candidate. In order to stabilize the fermionic singlet against a decay, we assume the dark sector particles to be charged under an unbroken stabilizing symmetry. We found five Dirac and five Majorana models divided into three different classes (labeled quarkphilic, leptophilic and amphiphilic) sharing similar phenomenological features. For each model, we studied DM relic density, spin-dependent and spin-independent direct detection in a hierarchical coupling set-up (Γ b Γ s ), as well as in democratic coupling set-up (Γ b ∼ Γ s ). Furthermore, we studied different mass-hierarchies between the DM and the other dark-sector particles within each model, ranging from coannihilation scenarios (κ = 1.01, 1.1) to non-coannihilation scenarios (κ = 5, 15).
Since these models typically offer a vast parameter space, we focused this study on the models' ability to solve the aforementioned anomalies in the major part of the parameter space under study. This is achieved by keeping the product of the second and third generation quark couplings Γ * s Γ b at the maximal value allowed by B-B-mixing, while Γ µ is left unconstrained.
Tables IV and V summarize the models' achievements in each set-up regarding their potential to provide simultaneous (Tbl. IV) or individual solutions (Tbl. V) to the R K and (g − 2) µ anomalies and the observed relic density under the premise of accordance with direct detection bounds.
We found that all Dirac DM models are excluded as viable solutions to the R K and (g − 2) µ anomalies by direct detection, as the DM-nucleon cross section receives a large contribution from the kinematically unsuppressed (v, v) currents in Z-boson exchange. This is a particularly strong statement, when we remind ourselves that the one-loop-type models addressing the R K anomaly usually require a large coupling to muons, which diminishes the relic density substantially, leading to significant relaxation of the direct detection bounds. Majorana DM models, on the other hand, can offer open parameter space due to the vanishing vector currents in these type of models.
Leptophilic Majorana models in principle allow for a solution to R K with underproduced DM for coannihilation scenarios in the region M DM 1000 GeV. The (g − 2) µ -anomaly on the other hand can only be solved for DM masses 290 GeV. An individual solution to DM only is generally possible in the low-Γ µregions of coannihilation scenarios. Note that this statement is even true for Dirac DM leptophilic models while still R K and (g − 2) µ cannot be explained 18 . We also have shown that a hierarchical flavor structure generally allows for larger DM masses and relaxes direct detection bounds in this context. RD   In contrast, quarkphilic and amphiphilic Majorana models conceptually allow for a simultaneous solution of both R K and DM. While our study only found a window for large DM masses M DM 1.5 TeV in a coannihilation scenario for the quarkphilic model bIIB Maj , the amphiphilic aIA Maj offers a low-M DM solution for a non-coannihilation set-up. Additionally, we found a small window in the TeV range for a coannihilation scenario, where ψ is a viable DM candidate without solving the R K anomaly 19 . Quarkphilic models also offer individual solutions to R K , (g − 2) µ and DM in democratic coannihilation scenarios. In the case of a hierarchical flavor structure, direct detection rules out all viable parameter space due to large, Γ b -driven contributions from quark loops to spin-independent DM-gluon scattering and spin-dependent DM-quark scattering. Amphiphilic models are not able to pose even an individual solution for (g − 2) µ but provide a number of solutions to the R K anomaly, including coannihilation and non-coannihilation scenarios.
As Tbl. IV suggests, the most interesting models studied in this work contain colored particles in the TeV-range, which raises the question of detectability at colliders. Collider searches in similar set-ups can exclude dark sector particle masses up to the TeV scale. These results, however, assume fixed values of DM couplings, which motivates a customized collider study. 19 Besides the low-massive solution, which we expect to be strongly constrained by a collider study, all remaining viable scenarios feature colored coannihilation. Those setups might receive sizable corrections to the relic density from non-perturbative effects such as the Sommerfeld enhancement and bound state formation arising in the annihilations of the colored particles, as discussed for instance [56] for a simplified colored coannihilations scenario. Those effects are not considered in MICROMEGAS and an inclusion lies beyond the scope of this work.  Note added: During the finishing stages of this work, the article [57] appeared that studies the problem discussed from a slightly different perspective. Both works discuss one-loop extension addressing B-anomalies in the context of t-channel DM. Out of the ten models discussed in this paper four are covered in the analysis of [57], namely aIA and bIIB for both Majorana and Dirac fermions, albeit with a more elaborate collider phenomenology.
Moreover, the DM analysis in both works is complementary, since the works study different parts of the parameter space and the fixed mass ratios of DM to the other dark sector particles employed here are particularly useful to study coannihilating scenarios while the approach in [57] is more suitable for collider phenomenology.
Finally, we point out here that the spin-dependent limits are especially relevant for the case of Majorana DM where they can be more constraining than the spin-independent searches. annihilation cross section, efficient conversions of DM, encoded in the ratio of equilibrium densities in Eq. (35), might alter the predictions of the relic density significantly in the case of small mass splittings.
More precisely, we address the interplay of the Yukawa coupling Γ i required to reproduce the observed relic density for a given mass ratio κ.
If we only consider the thermally averaged cross section of the direct DM annihilation in Eq. (37), we find that the cross section decreases for an increasing κ for both Majorana and Dirac DM. This corresponds to a decreasing cross section for an increasing mediator mass.
If we additionally consider the contribution of the ratio of equilibrium densities at the time of thermal freeze-out, x f ≈ 25, the κ dependence of the effective thermally averaged annihilation cross section results in σ eff v ∼ 1 (1 + κ 2 ) 2 1 + g non-DM g DM κ