Explaining Dark Matter and $B$ Decay Anomalies with an $L_\mu - L_\tau$ Model

We present a dark sector model based on gauging the $L_\mu - L_\tau$ symmetry that addresses anomalies in $b \rightarrow s \mu^+ \mu^-$ decays and that features a particle dark matter candidate. The dark matter particle candidate is a vector-like Dirac fermion coupled to the $Z^\prime$ gauge boson of the $L_{\mu}-L_{\tau}$ symmetry. We compute the dark matter thermal relic density, its pair-annihilation cross section, and the loop-suppressed dark matter-nucleon scattering cross section, and compare our predictions with current and future experimental results. We demonstrate that after taking into account bounds from $B_s$ meson oscillations, dark matter direct detection, and the CMB, the model is highly predictive: $B$ physics anomalies and a viable particle dark matter candidate, with a mass of $\sim (5-23)$~GeV, can be accommodated only in a tightly-constrained region of parameter space, with sharp predictions for future experimental tests. The viable region of parameter space expands if the dark matter is allowed to have $L_\mu-L_\tau$ charges that are smaller than those of the SM leptons.


Introduction
The presence of non-baryonic dark matter in the universe has been ascertained from a variety of observations corresponding to different time scales in the history of the universe. All such observations hinge, however, exclusively on gravitational effects associated with the dark matter, and no striking evidence for non-gravitational manifestations of dark matter is known to date. A broad experimental program is presently in place, aimed at discovering such non-gravitational effects with collider searches and with direct and indirect dark matter detection.
Collider searches for dark matter typically rely on so-called mono-X searches, i.e. missing energy plus some visible final state (see [1][2][3][4][5][6][7][8][9][10][11][12][13] for the latest analyses); such searches have the advantage of being sensitive to low dark matter masses, as opposed to e.g. current direct detection experiments which loose sensitivity for small masses below a few GeV. In addition to the usual mono-X studies, colliders might also give important information on the mediator that connects the dark and visible sector, for example in the case of vector mediators, as we will describe in more detail below.
In the class of models we investigate here, the dark matter is connected to Standard Model (SM) particles via a new gauge interaction that is also relevant to address putative new physics (NP) signals observed in rare B meson decays by the LHCb collaboration.
Specifically, we build on the model presented in Refs. [14,15], which is based on gauging L µ −L τ , the difference of leptonic muon number and tau number [16][17][18]. The model of [14,15] was proposed to explain an anomaly in the rare B → K µ + µ − decay mode observed by the LHCb collaboration in 2013 [19]. 1 The B → K µ + µ − anomaly persists in the latest experimental data [36,37], and is supported by measurements of the B s → φµ + µ − [38] decay and by the hint for lepton flavor universality violation in B + → K + + − [39]. Interestingly enough, the model in [14,15] necessarily predicts lepton flavor universality violation and is in excellent agreement with the latest model independent NP fits that take into account all relevant experimental data on rare B decays [40][41][42][43].
Beyond being able to explain the hints for lepton flavor universality violation in rare B decays, gauged L µ − L τ has a number of additional virtues. It has been considered as a solution to the anomaly in the g − 2 of the muon, and is used in models for neutrino masses. Unbroken L µ − L τ predicts one degenerate neutrino pair and a maximal atmospheric neutrino mixing angle, which is a good starting point for model building.
Extra matter could be charged under the L µ − L τ gauge symmetry. Here, we consider a minimal dark sector composed by the Z gauge boson associated to the L µ − L τ symmetry and by a vector-like Dirac fermion that is a dark matter candidate, and which we assume to be charged under the new local L µ − L τ symmetry. We investigate whether both the rare B decay anomalies and dark matter can be simultaneously accounted for, considering all relevant phenomenological constraints. The dark matter phenomenology is dictated by gauge interactions mediated by the massive new Z gauge boson. The mass of the Z can be generated either through a Stueckelberg mechanism, thereby leaving the L µ − L τ symmetry unbroken [44,45], or by a spontaneous symmetry breaking mechanism with the introduction of a scalar field, as described in Ref. [14]. We will assume that the scalar field, if it exists, is either sufficiently heavy to be integrated out and inessential to the phenomenology of the model, or that is has sufficiently suppressed interactions with SM particles, again leading to no observable consequences. In what follows we remain agnostic about the mechanism by which the Z gauge boson acquires mass, and we outline the region of parameter space which accommodates both dark matter and the rare B decay anomalies.
Notice that the construction under consideration here is primarily motivated by the notion of minimality in addressing both dark matter and B physics anomalies. As a result, while a full underlying model might reveal a connection between the ingredients of the model (e.g. the vector-like dark matter particle and the vector-like quarks that generate the effective couplings between the Z and the Standard Model quarks), none is postulated here. The construction of a complete model, including connections among neutrino masses, dark matter and vector-like quarks goes beyond the scope of this work.
Previous work has explored dark matter in models based on the L µ − L τ symmetry. Early studies of dark matter charged under L µ − L τ were motivated by the PAMELA positron fraction excess and the anomalous magnetic moment of the muon [46][47][48] and the galactic center excess [49]. Extended L µ − L τ setups have also been explored, containing right-handed neutrinos and additional light scalars [50] or additional gauge bosons [51]. More studies of L µ − L τ models with extended fermion and scalar sectors have been performed in [52,53]. A complementary study of the interplay of B physics anomalies and Dark Matter in models with gauged Baryon number can be found in [54].
Our work extends and complements previous studies in several notable directions: • We postulate a Dirac fermion as a dark matter candidate and discuss several bounds on the Z gauge boson stemming from neutrino trident production and collider searches in the context of a local L µ − L τ symmetry.
• We consider stringent limits stemming from indirect detection, specifically from CMB measurements on the energy injection at dark ages.
• We account for 1-loop dark matter-nucleon spin-independent scattering, and compare such predictions with the LUX and projected XENON1T experimental sensitivity. Despite the 1-loop suppression, we find remarkably strong constraints from direct detection experiments.
• Lastly, we outline the region of parameter space where both dark matter and the rare B decay anomalies can be simultaneously addressed in agreement with all existing bounds.
The paper is organized as follows: Sec. 2 introduces the model under consideration; Sec. 3 discusses general constraints on the Z parameter space, while the following Sec. 4 details on the connection with the B decay anomalies and Sec. 5 discusses the dark matter-related phenomenology. Finally, Sec. 6 describes and summarizes our findings, and Sec. 7 concludes.

The L µ − L τ Model
In the L µ − L τ model, the SM is augmented by a new Abelian symmetry, U (1) µ−τ , with a corresponding new massive gauge boson, Z . We remain agnostic as to how the Z acquires its mass, but we do assume that the physics connected to the mass generation is sufficiently decoupled and phenomenologically irrelevant for our discussion. The new boson couples to the SM fields via the covariant derivative D α = ∂ α + ig q µ−τ Z , with g the U (1) µ−τ coupling strength and q µ−τ the corresponding charge, according to the neutral current where q is a free parameter which quantifies the overall charge of the leptons under the L µ − L τ symmetry. i indicates the SM lepton doublets with flavor i. This interaction results into, Dark matter can be charged under this new gauge symmetry. A simple extension of the minimal model that incorporates a dark matter candidate is constructed by adding a vector-like Dirac fermion χ, singlet under the SM gauge group but charged under the new U (1) µ−τ symmetry. The fermion has a vector-like mass, m χ , and its coupling to the Z is given by where q χ is the dark matter charge under the U (1) µ−τ symmetry. The vector-like nature of the dark matter particle ensures the absence of triangle anomalies. In our analytic results we keep the dependence on the two charges q and q χ explicit. For our numerical results, however, we adopt q = 1, without loss of generality. The larger the q L µ − L τ charge, the smaller the g value needed to reproduce the same results. We will mainly concentrate on the case q χ = q , i.e. we assume that the Dirac fermion dark matter candidate features a "universal coupling" to the Z , equal to that of leptons. We will also comment in what follows on the potentially phenomenologically interesting case q χ q .

Constraints on the Z Parameter Space
A powerful probe of a Z vector boson based on gauging L µ − L τ is the process of neutrino trident production [55], i.e. the production of a µ + µ − pair in the scattering of a muon neutrino in the Coulomb field of a heavy nucleus. Integrating out the Z , the correction to the trident cross section can be concisely written as [14] where v = 246 GeV is the electroweak vacuum expectation value (vev) and s W ≡ sin θ W is the sine of the Weinberg angle θ W . Using the measurement of the trident cross section by the CCFR collaboration [56], one finds, for a given g gauge coupling, the following lower bound on the Z mass This bound is compatible with very light Z bosons as long as the coupling g is very small. For a light Z , with mass comparable to the neutrino momentum transfer in the trident reaction, the approximate expression in (3.1) breaks down. In this region of parameter space the Z has to be kept as dynamical degree of freedom. For m Z 10 GeV we use the results of a numerical evaluation of the trident cross section in the presence of a light Z from Ref. [55]. For such light Z the constraint in (3.3) is slightly relaxed. The region that is excluded by the trident measurements is shaded in red in Figs. 3 -5. We stress that this constraint is independent of the DM parameter space as defined by the parameters (q χ , m χ ).
An additional constraint on the parameter space of the L µ − L τ gauge boson comes from measurements of the Z → 4 branching ratio [14]. The ATLAS collaboration has measured [57] the fiducial branching ratio BR(Z → 4 ) = (4.2 ± 0.4) × 10 −6 , (3.4) using the full 7 and 8 TeV data set. This measurement is in good agreement with the SM prediction BR(Z → 4 ) SM = (4.37 ± 0.03) × 10 −6 . (3.5) An additional ATLAS Run I analysis [58] leads to very similar results. Also the CMS collaboration has performed a search for Z → 4 using 2.6 fb −1 of 13 TeV data [59]. The 13 TeV result, however, has a still larger (statistical) uncertainty, if compared to the ATLAS search based on the full 8 TeV data set. We do not attempt a statistical combination of these results, but will use the measurement (3.4) in the following. Our Z contributes to the process with the 4 muon final state. Despite the fact that the ATLAS search has not been optimized to specifically constrain the Z scenario, interesting bounds on the gauge coupling g can be set for 4 m Z /GeV 70 where the three-body decay Z → µ + µ − Z is open. For larger Z masses the phase space starts to close, while for lower Z masses, the experimental acceptance becomes too small: the ATLAS analysis requires, in fact, two independent pairs of leptons with invariant masses above 5 GeV, implying that in order to get a non zero acceptance for m Z < 5 GeV, one would need to "mispair" the leptons, since, otherwise, one pair would have an invariant mass m µµ = m Z < 5 GeV. Secondly, going to lower and lower Z masses, the two leptons coming from the Z decay become more and more collimated, failing more often the isolation criterium employed in [57].
The region in the m Z -g plane that has been probed by the Z → 4 measurement is shown in gray, in Figs. 3 -5. The most stringent limit on the gauge coupling, g 0.015, can be set at around the mass m Z 10 GeV. It might be possible to substantially extend the collider reach in the coming years of the LHC using targeted new searches for Z → 2µZ , Z → 2µ and for Z → 2µZ , Z → νν [60]. Notice that the LHC bound on the Z depends slightly on the DM mass: if the Z is heavier than two times the DM, then the Z will also decay to a DM pair, weakening in such a way the bound coming from the measurement of Z → 4 . We take into account this effect in our numerical results. Finally, we point out that the BaBar collaboration has recently performed a search for muonic forces measuring the cross section for the process with the Z decaying into µ + µ − [61]. The search constrains regions of Z parameter space with 2m µ < m Z 8 GeV. The corresponding constraints on g are slightly stronger than the trident constraints for m Z 4 GeV. Such low Z masses are outside of our main region of interest. Therefore the resulting bound will not be presented in what follows.

Rare B Decay Anomalies
Over the last few years, various anomalies in rare B meson decays have been reported by the LHCb Collaboration [62]. Such anomalies include discrepancies in the B → K * µ + µ − angular distribution [19,36] (confirmed by a recent Belle result [37]), a reduced branching ratio for the decay mode B s → φµ + µ − [38], as well as a hint for lepton flavor universality violation in B → Kµ + µ − vs. B → Ke + e − [39]. Assuming that hadronic uncertainties in the corresponding SM predictions are estimated in a sufficiently conservative way, global fits of the combined rare B meson decay data show a strong preference for new physics in b → sµµ transitions, while b → see transitions appear compatible to the corresponding SM predictions [40][41][42][43]. The best description of the data is obtained by new physics in the form of a four fermion contact interaction with C NP 9 = −1.07 and V is the CKM matrix. The Z arising from the L µ − L τ gauge symmetry is ideally suited to address the B physics anomalies as it has the required vector couplings to muons and does not couple, by construction, to electrons. As shown in Ref. [14], introducing effective flavor-changing couplings of the Z gauge boson to left handed quarks and integrating out the Z leads precisely to the contact interaction in Eq. (4.1).
The Wilson coefficient C NP 9 is determined by three parameters: the Z mass, m Z , the Z coupling to muons, q g , and its flavor violating b ↔ s coupling. In the following we choose m Z and g as free parameters and set the flavor-violating coupling such that the rare B decay anomalies are explained, i.e. such that the best fit value for C NP 9 = −1.07 from [41] is reproduced. This is possible as long as the Z boson is sufficiently heavy compared to the B mesons. A m Z below the mass of B mesons would lead, in fact, to a resonance in the dimuon invariant mass spectrum in rare B decays which is strongly constrained [63]. To ensure a consistent explanation of the B anomalies we will demand the conservative bound m Z 10 GeV. Note that for such Z masses, the effective operator approach of eq. (4.1) is fully justified. The region of di-muon invariant mass squared that is most relevant for the rare B decay anomalies is 1 GeV 2 < q 2 < 6 GeV 2 . For m Z 10 GeV, corrections are therefore at the few percent level at most.
Any explanation of the B decay anomalies based on Z bosons is subject to additional constraints from B s meson oscillations. Indeed, integrating out the Z does not only lead to contributions to the b → sµ + µ − decays. Once the flavor violating b ↔ s coupling is fixed to explain the anomalies, it necessarily leads also to corrections to B s meson oscillations. We find the following modification of the B s −B s mixing amplitude M 12 with the SM loop function S 0 2.3 [64]. Allowing for |M Z 12 |/|M SM 12 | 15% [65, 66] and setting C NP 9 = −1.07 [41], as favored by the anomalies in B meson decays, we obtain the following upper bound on the Z mass This bound in shown in Figs. 3 -5 in magenta. In summary, the parameter space that allows to address the B decay anomalies, while being consistent with the trident and Z → 4 constraints, is characterized by a lower bound on the ratio m Z /(q g ) coming from neutrino tridents, an upper bound on m Z /(q g ) coming from B s meson oscillations and a lower bound on m Z from the dimuon invariant mass distribution in rare B decays 540 GeV m Z /(q g ) 4.9 TeV , m Z 10 GeV . (4.4) Outside this window of favored parameter space, the model cannot address the B decay anomalies. We conclude this section noting that this parameter region is completely independent of the DM parameter space (q χ , m χ ).

Dark Matter Phenomenology
In this section we discuss the phenomenology of the dark matter in our model. It is important to first notice that the Z gauge boson arising from the L µ − L τ local symmetry dictates a dark matter phenomenology rather different from other realizations of a "dark Z portal" (see e.g. [67][68][69][70][71][72][73][74][75][76][77][78][79][80][81][82][83][84]). For example, the dark matter pair annihilation exclusively yields muon, tau, and neutrino pairs, as well as Z boson pairs, if kinematically allowed (i.e. if m χ > m Z ). We show the schematic relevant Feynman diagrams in Fig. 1. This feature affects both the relic density and the indirect detection properties of the model. Secondly, scattering off of nuclei (direct detection) occurs only via loop induced dark matter-nucleon interactions (see Fig. 2). Most of the previous results obtained in the literature are therefore not directly applicable here.

Relic Density
In the early universe, the dark matter particle was in thermal equilibrium with SM particles through the interactions shown in Fig. 1. As the universe expanded and cooled down, the expansion rate eventually became equal to the interaction rate, leading to the freeze-out of the dark matter particles, and to a relic dark matter population. Since in our model the interactions of the dark matter particles with SM particles are dictated by s-channel (velocity-independent) processes, away from the m χ m Z /2 resonance, the abundance is directly tied to the annihilation times relative velocity of the annihilating DM particles today. In the case of charged leptons in the final state, such pair-annihilation cross section times relative velocity is given by In the case of neutrinos in the final state the expression is very similar for each neutrino flavor ν = ν µ , ν τ . Finally, in the case of final state gauge bosons Z , relevant when m χ > m Z . The quasi-on-shell Z exchange for m χ m Z /2 gives rise to a resonant enhancement of the cross section. In this region of parameter space, Eqs. (5.1) and (5.2) have to be changed taking into account the Z width. In our numerical analysis this is taken into account using the micrOMEGAs code [85,86]. From Eqs. (5.1) -(5.3) we learn that annihilation into neutrinos and charged leptons occurs at similar rates, and is also comparable to the annihilation rate into Z gauge boson pairs if m χ > m Z . Given the annihilation cross section above, it is straightforward to obtain the dark matter thermal relic density, which is approximately with the Planck mass M pl = 1.22 × 10 19 GeV, x F the inverse temperature at freezeout in units of the dark matter particle mass, and g the number of relativistic degrees of freedom at freeze-out (around 90 for freeze out temperatures of (5-80) GeV [87,88]). Notice that for the actual numerical evaluation of the dark matter thermal relic abundance we use the micrOMEGAs code [85,86]. In Figs. 3 -5 we show in green the curves that reproduce the observed thermal relic abundance Ω χ h 2 0.12 [89] as a function of the g coupling and the Z mass, for a given value of the particle dark matter mass. Above these curves our DM candidate is under-abundant. The figures illustrate that there are three distinct regimes: (i) for m χ m Z , the relevant annihilation cross section scales as the combination g 4 /m 2 χ (for given choices of q χ , q f ): the correct relic density thus singles out one value of g for a given dark matter particle mass m χ ; (ii) for m χ m Z /2 the resonant regime sets in, driving the g producing the right relic density to suppressed values; (iii) for m χ m Z the cross section scales as the combination g 4 m 2 χ /m 4 Z , and thus the right relic density selects a value for the ratio g /m Z , for a given value of m χ .
The regions of correct thermal relic density on the plane defined by the g coupling versus the dark matter mass for a given Z mass are shown in green in Fig. 6. The main features of the curves can be again simply understood from the approximate analytic form for the pair-annihilation cross sections in Eqs. (5.1) -(5.3). The "funnel" at low values of g corresponds to the resonant annihilation mode, while the large g vertical asymptote to the m χ m Z regime described above and, finally, the constant g 2 /m χ region to the asymptotic m χ m Z regime (which onsets fully at larger m χ than those in the figure).

Indirect Detection
Several probes exist of the late-time, low relative-velocity dark matter pair-annihilation cross section, generically known as "indirect dark matter detection". The characteristic feature of the model under consideration is the production of abundant hard leptons, as illustrated by the annihilation diagrams shown in Fig. 1. As such, one possible constraint stems from cosmic-ray and gamma-ray data. In particular, the AMS-02 measurement of the cosmic ray positron fraction can be used to set constraints on dark matter annihilating to leptons [90][91][92][93][94]. Dark matter below ∼ 100 GeV annihilating into muons is in principle subject to such constraints. However, utilizing charged cosmic rays to set constraints on dark matter annihilation inevitably involves significant uncertainties from propagation and energy losses in the Galaxy, and and, in the case of light dark matter, m χ 15 GeV, also from solar modulation. Though gamma-rays offer a promising search channel for dark matter annihilation, in the case of leptonic interactions, they are not the best probe. A more model-independent and robust way to probe the pair-annihilation rate is to use the effects of energy injection from dark matter annihilation on the cosmic microwave background (CMB). Specifically, dark matter annihilation at redshifts z ∼ 1000 results in energy injection which heats and ionizes the photon-baryon plasma, significantly perturbing the ionization history in a way that can be constrained by measurements of the CMB temperature and polarization angular power spectra [95]. In our model, where the SM annihilation final states largely involve charged leptons, a sizable fraction, and in some case almost all of the energy is deposited in energetic electrons and positrons, the species with the highest effective deposited power fraction f ef f [95].
(5.6) If the Z boson is sufficiently light also χχ → Z Z annihilations are possible. Neglecting again phase space effects and setting q χ = q , we find approximately with the Z bosons decaying back to neutrinos, muons and taus, if kinematically allowed.
In the case of 100% annihilation into muons, for which f ef f 0.2, the latest results from Planck [89] solidly exclude the canonical annihilation cross section of 3 × 10 −26 cm 3 /s for dark matter masses below ∼ 20 GeV [89], when no velocitydependence exists in the pair-annihilation cross section of the Dirac dark matter. Pair-annihilation to taus, which also includes hadronic final states, leads to similar results since it has only a slightly lower f ef f . Annihilation to four muons (which is relevant in the case of χχ → Z Z → 4µ) has also been shown to correspond to f ef f ∼ 0.2 [103], leading therefore to similar conclusions. Annihilations into neutrinos produce negligible effects to the CMB power spectrum. In the latter scenario, neutrino [104,105] and gamma-ray detectors provide a better probe [106], but these are still far from the canonical cross section, thus placing no relevant constrain to our model.
The portion of the pair-annihilation cross section giving rise to the right amount of thermal relic dark matter (green curve) which is excluded by current Plack data is shaded in blue and labeled by CMB in Fig. 3. In addition to the current Planck results, we will also quote a forecast for a cosmic-variance-limited experiment with similar angular resolution which might improve the current sensitivity by a factor of four [89]. The corresponding region is shaded in dark blue and labeled as CMB proj. in Figs. 3 -6. The canonical annihilation cross section is excluded for m χ 20 GeV for 100% annihilation into muons using current CMB data. However, in our case, a sizeable fraction of the annihilation rate goes into neutrinos that yield negligible constraints. For this reason we can see that for m χ = (10 − 15) GeV, Figs. 3-5, only the projected CMB sensitivity applies. It is important to emphasize that, at the resonance, the direct relation between the annihilation cross section today and at the time of freezeout might fail, as pointed out in Ref. [107,108], due to thermal effects in the early universe, explaining why the resonance region in the figures is not excluded.

Direct Detection
Since the DM particle does not directly couple to quarks, scattering off of nuclei occurs only through a loop of charged leptons that couple to photons which in turn couple to protons, as illustrated in Fig. 2 2 . The WIMP-nucleon scattering cross section is thus proportional to the nucleus electric charge. Moreover, since the Z boson has vector-like interactions with the charged leptons and the dark matter, the scattering cross section is not velocity dependent. The left diagram in Fig. 2 is the most relevant, while the other two are 2-loop suppressed and therefore negligible. In our numerical calculations we include only the former. Adapting the results obtained in [109] for leptophilic dark matter to our model, we find the WIMP-nucleon scattering cross section to be where µ N = m N m χ /(m N + m χ ) is the WIMP-nucleus reduced mass, m N the nucleus mass, and Z and A the atomic and mass numbers, respectively. Since we will compare our theoretical predictions with current limits from Xe-based experiments, Z = 54 and m N 129 GeV. The dependence on the logarithm of the lepton masses squared can be understood as a leading log approximation to the RGE induced kinetic mixing between the Z and the photon in the running from the tau mass down to the muon mass.
The most stringent limits on the WIMP-nucleon spin-independent scattering cross section stem at present from the LUX experiment [110], which exclude σ SI 2× 10 −10 pb for m χ ∼ 50 GeV, surpassing earlier results [111,112] after achieving better calibration, event-reconstruction, background modeling, and more live-days. Similar sensitivity has been achieved by PANDAX-II [113]. We present the current LUX and the projected Xenon1T sensitivities in solid green in Figs. 3-5, for a variety of choices for the dark matter particle mass, set to 5, 10 (Fig. 3), 15, 50 and 100 GeV (Fig. 4).
vector-like quarks, which mix with SM ones, inducing a small flavor violating Z bs coupling. While a flavor-violating coupling cannot lead to any appreciable contribution to the direct detection cross section, it might be accompanied by flavor diagonal couplings to quarks. The size of the flavor diagonal couplings is model dependent and we will not consider their effects here.
The figures also show the projected sensitivity from the XENON1T experiment with 2 years of data [114], expected to improve the current LUX limits by two orders of magnitude (dashed green line). To obtain these curves we assume canonical values for the local dark matter density and velocity distribution throughout the plots.
The figures illustrate that direct detection provides remarkably strong limits in spite of being loop-suppressed in our model. Direct detection constraints exclude dark matter masses above ∼ 15 GeV, allowing only the peak corresponding to the Z resonance. For m χ < 15 GeV, direct detection limits weaken, because for lower masses the recoil energy is sufficiently small and close to the energy threshold of the experiments, degrading the corresponding sensitivity. This conclusion depends on the fact that we have fixed q = q χ . In the next section we will also depart from this equality and address the consequences. The corresponding results are collected in Fig. 5.

The Global Picture: Dark Matter and B Anomalies
In this section we combine all findings described in the previous sections and outline the region of parameter space which can simultaneously accommodate thermal relic dark matter and the rare B decay anomalies. The region in the (m Z , g ) parameter space that accommodates the B physics anomalies corresponds to the central diagonal white regions in Figs. 3-5. The top red and bottom magenta shaded regions are ruled out by neutrino trident production and B s mixing, respectively. The gray contour delimits the region ruled out by the measurement of the Z → 4µ decay width. Note that the region favored by the flavor anomalies ends at ∼ 10 GeV, as discussed in Sec. 4. The light green curves indicate (m Z , g ) combinations that reproduce the right relic abundance (Ω χ h 2 0.12). Relic abundance curves with the light (dark) blue contours are excluded by limits from current (projected) CMB data. Diagonal solid (dashed) green lines indicate the LUX (the projected XENON1T) bounds: the region to the top left of those lines is excluded.
In Fig. 3 we show results for m χ = 5, 10 GeV and charge q χ = 1. For m χ = 5 GeV, current and projected direct detection limits (green curves on the top left of the plots) are weak due to the experimental energy threshold. The correct relic abundance can be achieved close to the resonance region, without being excluded by the constraints from the neutrino trident production. The CMB limit probes a sizable part of parameter space that yields the right abundance. Forecasted CMB data constrain the model a bit further. In the m χ = 5 GeV case we are left with a small region of parameter space with a Z mass of m Z 10 GeV and a L µ − L τ gauge coupling of g ∼ 10 −3 − 10 −2 .
For m χ = 10 GeV, only the projected CMB sensitivity applies. Direct detection limits from LUX and projected limits from XENON1T probe substantial parts of the parameter space of the model where one can accommodate both the LHCb anomaly  . The charge of the DM is always fixed to q χ = 1. The region able to explain the rare B decay anomalies, while being not ruled out by other experiments, is the central white diagonal band, for m Z 10 GeV. The light green curve corresponds to the parameters producing the right thermal relic dark matter abundance, with the portions shaded in blue excluded by CMB data at present (light blue) or in the future (darker blue). The shaded magenta and red regions are excluded by B s mixing and neutrino trident production observables, respectively. The gray contour delimits the region ruled out by the measurement of the Z → 4µ decay width. The dark green solid diagonal line represents the LUX-2016 direct dark matter detection bound, with the parameter space to the upper left being excluded. The dashed dark green diagonal line corresponds to the expected XENON1T sensitivity with 2 year data. and the relic density. Notice that CMB limits weaken at the resonance because, at the resonance, sufficiently small couplings can still reproduce the right abundance, however such small coupling lead to small annihilation cross sections today. This   trend can be seen for m χ = (5 − 15) GeV (see the upper panel of Fig. 4 for the latter case). Fig. 4 focuses on dark matter masses m χ = 15, 50, 100 GeV. Notice that direct detection limits (green curves on the top left of the plots) are significantly stronger than in the 5 GeV case, becoming most relevant for a 50 GeV dark matter mass; XENON1T projected sensitivity can entirely probe the LHCb favored region for m χ = 15 GeV. The entire favored parameter space for m χ = 50 GeV and for m χ = 100 GeV is already excluded by current LUX bounds combined with B s mixing.
We find that the current LUX bounds leaves viable parameter space only for m χ 23 GeV.
In summary, in the regime with q l q χ , the dark matter mass window of m χ (5 − 23) GeV is favored. For smaller dark matter masses the Z required to accommodate the right relic abundance is too light to explain the B decay anomalies. For larger dark matter masses direct detection and B s mixing exclude the entire favored parameter space. Upcoming XENON1T data might significantly test the model, leaving only masses around m χ (5 − 10) GeV open. The latter can be later probed with further improvements from B s mixing constraints, and the LZ direct detection experiment [115].
These conclusions change if we depart from the assumption that q l q χ . In Fig. 5 we show the results for q χ = 1/6, keeping q l = 1. The indirect detection limits are essentially the same since the smaller q χ the larger the gauge coupling needed to obtain the right relic abundance, yielding no impact in the overall annihilation cross section needed to get the relic density. On the order hand, direct detection limits are weaker ameliorating the scenario specially for heavy masses. Indeed, for the examples m χ = 500 GeV and 1 TeV, one can see viable regions of parameter space in the center and bottom panels of Fig. 5 that are able to simultaneously accommodate dark matter and the LHCb anomaly.
Our final Fig. 6 provides an orthogonal view of the parameter space in the (m χ , g ) plane, with the mass of the Z set to m Z = g × 4500 GeV which is inside the LHCb favored region for both choices of q χ : q χ = 1 (upper panel) and q χ = 1/6 (lower panel). In this regime, CMB limits do not touch the relic density curves (green) but do probe the funnel region with larger annihilation cross sections. In the figures we show these limits in dark blue. For q χ = 1/6, the funnel region is rather narrow and for this reason the CMB projected limit is simply a line. The vertical red lines show contours of constant spin-independent dark matter-nucleon cross section σ SI . Notice that since σ SI ∝ g 4 /m 4 Z and the ratio g /m Z is kept fixed in the figure, the cross section only depends on the dark matter mass through the reduced mass, µ N (see Eq. (5.8)). Having in mind that current LUX limit excludes σ SI 4 × 10 −9 pb for m χ = 10 GeV and σ SI 2.5 × 10 −10 pb for m χ = 100 GeV, we conclude that only light dark matter can circumvent all existing constraints and accommodate the B decay anomalies if q χ = 1, q = 1 (upper panel of Fig. 6). For the Z mass of m Z = g × 4500 GeV we find the bound m χ 22.9 GeV as indicated by the dashed vertical line in the plot.
In the bottom panel, we set q χ = 1/6, q = 1, weakening direct detection limits. In this case, the entire parameter space is consistent with the data and can address the B decay anomalies.

Conclusions
In this paper we proposed a new physics setup with dark matter charged under a L µ − L τ gauge symmetry. Specifically, we studied if one could address at the same time anomalies observed in rare B meson decays by the LHCb collaboration and have a viable thermal dark matter candidate. Our model hinges upon the (automatically anomaly-free) L µ − L τ local gauge symmetry, U (1) µ−τ , with a new corresponding massive gauge boson Z . Dark matter is a vector-like Dirac fermion charged under the L µ −L τ symmetry and neutral under the SM gauge symmetries. This setup leads to an unusual and novel dark matter phenomenology, which we described in detail.
The dark matter relic abundance is mainly set by annihilation into muon, tau, and neutrino pairs, through s-channel Z exchange. We studied indirect detection limits stemming from distortions of the CMB power spectrum, and direct detection using current LUX2016 and XENON1T projected sensitivities. Despite the fact that dark matter nucleus scattering is loop-suppressed, we found that direct detection experiments do put remarkably strong constraints on the model parameter space. The correct relic abundance can only be obtained close to the resonance region m Z 2m χ . Combining the direct detection constraints with constraints from B s meson oscillations leaves only a restricted window of parameter space where both B decay anomalies and the correct relic abundance can be explained: the dark matter mass is 5 GeV m χ 23 GeV; the Z mass is m Z 2m χ ; the gauge coupling is in the range 2 × 10 −3 g 10 −2 . Substantial parts of this parameter space can be probed by the expected sensitivities of future direct detection experiments like XENON1T and LZ.
The viable parameter space can be extended significantly by allowing the dark matter to have L µ − L τ charges that are smaller than the charges of the SM leptons. In this case direct detection constraints are weakened, and their upper bound on the dark matter mass disappears. A correct relic abundance still requires a Z mass that is close to the resonance region. Prospects for probing such a scenario with future direct detection experiments are excellent.

References
[1] ATLAS Collaboration, M. Aaboud et. al., Search for dark matter produced in association with a hadronically decaying vector boson in pp collisions at √ s=13 TeV with the ATLAS detector, 1608.02372. [7] CMS Collaboration, Search for dark matter and graviton produced in association with a photon in pp collisions at √ s = 13 TeV, 2016. CMS-PAS-EXO-16-039.
[9] CMS Collaboration, Search for dark matter in final states with an energetic jet, or a hadronically decaying W or Z boson using 12.9 fb −1 of data at √ s = 13 TeV, 2016. CMS-PAS-EXO-16-037.
[10] CMS Collaboration, Search for dark matter in association with a Higgs boson decaying into a pair of bottom quarks at √ s = 13 TeV with the CMS detector, 2016. CMS-PAS-EXO-16-012.