Flavourful $Z'$ portal for vector-like neutrino Dark Matter and $R_{K^{(*)}}$

We discuss a flavourful $Z'$ portal model with a coupling to fourth-family singlet Dirac neutrino dark matter. In the absence of mixing, the $Z'$ is fermiophobic, having no couplings to the three chiral families, but does couple to a fourth vector-like family. Due to mixing effects, the $Z'$ gets induced couplings to second family left-handed lepton doublets and third family left-handed quark doublets. This model can simultaneously account for the measured $B$-decay ratios $R_{K}$ and $R_{K^*}$ and for the observed relic abundance of dark matter. We identify the parameter space where this explanation is consistent with existing experimental constraints from dark matter direct and indirect detection, LHC searches, and precision measurements of flavour mixing and neutrino processes.


Introduction
Recently, the phenomenological motivation for considering non-universal Z models has increased due to mounting evidence for semi-leptonic B decays whose rates and differential distributions are inconsistent with those predicted by the Standard Model (SM) [1][2][3]. In particular, the LHCb Collaboration has reported a number of deviations from µ-e universality in B → K ( * ) l + l − decays. The ratios of µ + µ − to e + e − final states: R K [4] and R K * [5] are observed to be about 70% of their expected values, each displaying a 2.5σ deviation from the SM. Combining that with the input from other b → s + − processes, the SM is disfavored by 4 to 5 standard deviations [6,7].
The R K and R K * anomalies could be the first evidence of new physics. A number of recent phenomenological analyses, see e.g. [6][7][8][9][10][11][12][13][14], conclude that these data can be well fit when the low-energy Lagrangian below the weak scale contains a new physics operator of the C NP 9µ = −C NP 10µ form: (1.1) In a flavourful Z model, the new physics operator in eq. (1.1) will arise from tree-level Z exchange: G bsµ = − g bs gµµ M 2 Z , where g bs is the flavour-violating Z coupling to left-handed b-and s-quarks, and g µµ is the couplings to left-handed muons. There is already a vast literature discussing the Z explanation of the B-anomalies and phenomenological constraints on the parameter space of such models, see e.g. . In realistic models of this kind, the coupling g bs is strongly constrained by precision measurements of the B s meson mass difference. Taking that into account, one can derive the constraint M Z 1.2g µµ TeV, implying that M Z must be close to the weak scale in weakly coupled models. The corollary is that the Z is in the correct mass range to act as mediator between the SM and thermally produced dark matter [46][47][48][49][50][51][52]. In this paper we further pursue this direction, and discuss a Z model that can account for the B-anomalies and, simultaneously, explain the observed relic abundance via a weakly interacting massive particle (WIMP) communicating with the SM through the same Z .

JHEP08(2018)061
We follow ref. [53], which introduces a fourth vector-like family with non-universal gauged U(1) charges. The idea is that the Z couples universally to the three chiral families, which then mix with the non-universal fourth family to induce effective nonuniversal couplings in the physical light mixed quarks and leptons. Such a mechanism has wide applicability, for example it was recently discussed in the context of F-theory models with non-universal gauginos [54]. Two explicit examples were discussed in [53]. Firstly an SO(10) → SU(5) × U(1) X model, where we identified U(1) ≡ U(1) X , which however was subsequently shown to be not consistent with both explaining R K * and respecting the B s mass difference [55]. Ref. [53] also discussed a fermiophobic model where the gauged U(1) charges are not carried by the three chiral families, only by fourth vector-like family. In the absence of mixing, the Z is fermiophobic, having no couplings to the three chiral families, but does couple to a fourth vector-like family. Due to mixing effects, we shall suppose that the Z gets induced couplings to second family left-handed lepton doublets (containing the left-handed muon and its neutrino) and third family left-handed quark doublets (containing the left-handed top and bottom quarks). Including only such couplings is enough to address the B-anomalies, in analogy to related scenarios where new vector-like fermions mix with the SM ones [20,23,26,28,31,38,40,45]. In addition, this set-up provides a natural WIMP dark matter candidate: the neutrino residing in the fourth family. We are interested in the parameter space of this model where both B-anomalies and the relic abundance of dark matter are simultaneously explained. We show that this can be achieved without conflicting a myriad of direct and indirect dark matter constraints as well as experimental constraints such as B s mixing, LHC searches, neutrino trident, and so on. The requirement to satisfy all these constraints in a natural way points to a specific corner of the parameter space, with 300 GeV m Z 1 TeV, dark matter heavier than a TeV, and a narrow range of possible Z couplings. This paper is organized as follows. In section 2 we define our gauged U(1) model with a vector-like fourth family. The Z couplings relevant for the subsequent analysis are summarized in eq. (2.28). In section 3 we discuss the constraints these parameters need to satisfy in order to address the B-anomalies without conflicting other experimental results. In section 4 we turn to the dark matter sector, and identify the masses and couplings of the vector-like fourth family singlet Dirac neutrino which lead to a correct relic density, while evading all indirect and direct searches so far. Our main results are contained in section 5, where we put together the requirements imposed by the B-anomalies and by the relic density, and identify the viable parameter space where both are satisfied.

The model
We consider a model in which, in addition to the SM with the usual three chiral families of left-handed quarks and leptons, including the right-handed neutrinos, we add a dark U(1) gauge symmetry and a fourth vector-like family of fermions. The idea is to have the SM quarks and leptons neutral under the U(1) while the vector-like family has the SM quantum numbers and is charged under the U(1) , leading to a dark matter candidate and flavourchanging Z operators after the vector-like fermion mass term mix with the SM fermions.

JHEP08(2018)061
Field Representation/charge  Table 1 shows all the particle content and their corresponding representations and charges. The non-universal U(1) charges forbid mixing between the fourth family and the chiral families via the usual Higgs Yukawa couplings. Therefore, we need to add new singlet scalars, with appropriate U(1) charges, to generate mass mixing of quarks and leptons with the vector-like family. The U(1) is broken by the VEVs of the new Higgs singlets φ ψ to yield a massive Z .
The Higgs Yukawa couplings of the first three chiral families can be written in a 4 × 4 matrix notation L Yukawa = y uQ LH u R + y dQ L Hd R + y eL L He R + y νL LH ν R + h.c. , (2.1) whereH = iσ 2 H * and y u , y d , y e , y ν are 4 × 4 matrices with the fourth row and columns consisting of all zeros, since the fourth family does not couple to the Higgs doublets. The U(1) charges allow Yukawa couplings between the singlet fields φ, the fourth familyψ 4 and the first three chiral families ψ i . Furthermore, there is an explicit mass term between the opposite chirality fourth family fields ψ 4 andψ 4 , where i = 1, . . . , 3.

JHEP08(2018)061
The fourth-family vector-like singlet neutrinos ν R4 ,ν L4 are special since we don't have a singlet field φ ν that couples them to the other families, which is why such terms are absent in the above equation. This implies that ν R4 ,ν L4 are absolutely stable, with their stability guaranteed by an unbroken global U(1) ν R4 and, since they do not carry any Standard Model quantum numbers, they may play the role of dark matter. Note that we also impose lepton number conservation U(1) L for all four families of leptons which forbids Majorana mass terms. Hence all neutrinos (including those in the fourth vector-like family) will have purely Dirac masses. 1 After the singlet scalar fields φ obtain a non-zero vacuum expectation value (VEV), we may rewrite the Lagrangian in terms of new mass parameters M Q i = x Q i φ Q , similarly for the other mass parameters, such that where α = 1, . . . , 4. We may diagonalize the mass matrix before electroweak symmetry breaking, when only the fourth family is massive The prime states for the heavy mass basis where only the fourth family has explicit vectorlike Dirac mass terms and it's related to the original charge basis by unitary mixing matrices, while for the neutrino statesν L4 and ν R4 the original and the mass basis coincides. In this basis, the Yukawa couplings in eq. (2.1) become where This shows that there is a coupling between the heavy fourth family and the Higgs due to their mixing with the first three chiral families. However, this coupling will be small since the original y u , y d , y e , y ν contain zeroes in the fourth row and column and they are mixing suppressed. Therefore, we can integrate out the fourth family and look at the low energy effective theory by simply removing the fourth rows and columns of the primed Yukawa matrices in eq. (2.6). The three massless families, below the heavy mass scale, are described by

JHEP08(2018)061
where and i, j = 1, . . . , 3. The Yukawa matrices for the quarks and charged leptons can be now diagonalized (2.10) The unitary CKM matrix is then given by In the case of neutrinos, since we are forbidding Majorana masses, the light physical neutrinos have Dirac mass eigenvalues given by, The lepton mixing matrix or PMNS matrix can be constructed from the transformations in eqs. (2.10) and (2.12) To look at the Lagrangian involving the SM gauge couplings, we emphasize that all the four families have the same charges under the SM. The unitary transformations in eq. (2.5) cancel as in the usual GIM mechanism and the gauge couplings in the heavy mass basis remains the same as in the SM. After integrating out the fourth family and electroweak symmetry is broken, and the light Yukawa matrices are diagonalised, the couplings to the W ± gauge bosons are (2.14) where g 2 is the usual SU(2) L gauge coupling. For the couplings to the Z gauge boson, the same happens, the charges are the same for the fourth families and the transformations in eq. (2.5) cancel, such that in the heavy mass basis, after electroweak symmetry breaking, we are left with and

JHEP08(2018)061
The electric charge of the fermions is denoted by Q and t 3 are the eigenvalues of σ 3 /2. The couplings to the Z boson are flavour diagonal, even after diagonalization of the light fermion mass matrices, due to the unitary transformations cancelling. The interactions will be the same as in eq. (2.15), replacing the fields ψ α by their three family mass eigenstates.
In the case of the couplings to the Z gauge bosons, we have non-universal couplings that lead to flavour changing. In the original basis, after the U(1) symmetry is broken, we have diagonal gauge couplings between the massive Z gauge boson and the four families where, In addition there are the fourth family couplings involving the opposite chirality states ψ 4 . Using the transformations in eq. (2.5), we get the Z couplings in the diagonal heavy mass basis and similarly with Q → L, etc. Ignoring phases, these matrices can be parametrized as where s ij and c ij refer to sin θ ij and cos θ ij (we have also suppressed the superscript in the angles s Q 14 → s 14 for simplicity). Since the U(1) charges differ for the fourth family, the unitary transformations do not cancel and the matrices D Q , etc., are not generally diagonal. Therefore, Z exchange can couple to light families of different flavour.
We are interested in thesbZ andμµZ couplings, needed for the R K anomaly. Assuming that only the mixing angles θ Q L 34 and θ L L 24 are different from zero 2 the mixing mass matrices become 2 A more natural possibility would be to assume that the new vector-like fermions have a large mixing only with the 3rd generation of the SM doublet, that is with taus instead of muons. Then the coupling to muons could arise due to a mixing between the SM charged leptons, as in [53]. However, explaining the B-meson anomalies in such a set-up runs in conflict with the strong bounds from non-observation of τ → 3µ.

JHEP08(2018)061
while the rest of them being zero. In the low energy effective theory, after integrating out the fourth heavy family, the Z couplings to the three massless families of quarks and leptons are where Q L3 = (t L , b L ) and L L2 = (ν µL , µ L ). Using now the diagonalization of the Yukawa matrices in eq. (2.10), we can expand the primed fields in terms of the mass eigenstates, For simplicity, we assume that the charged lepton mass matrix is diagonal so that we may drop the primes on the muon field so that µ L = µ L . Under this assumption, in the lepton sector, the Z only couples to muon mass eigenstates µ L and muon neutrinos ν µL , where the latter are related to neutrino mass eigenstates by the PMNS matrix, Given the hierarchies of the CKM matrix, we will assume similar hierarchies of the rotation matrix elements: The vector-like neutrino ν 4 is not charged under the SM and it is considered as a dark matter candidate. The portal that allows it to annihilate into ordinary matter is the Z mediator. The explicit coupling between the Z and the dark matter candidate ν 4 is where the Dirac dark matter field is given by We finish this section by summarizing all non-SM interactions that will later be relevant for our phenomenological analysis, introducing the notation that we shall subsequently use: where |V ts | ≈ 0.04 is the 3-2 entry of the CKM matrix, as otherwise unnatural cancellations would be required. It follows that |g bs | |V ts g bb |; in the following for simplicity we assume g bs = V ts g bb , and that g bb and g µµ have the same sign. Thus, the relevant parameter space is 5-dimensional: 3 couplings (g bb , JHEP08(2018)061 g µµ , g νν ) and 2 masses (M Z and the dark matter mass m ν ). From the theory point of view these are all essentially free parameters, although one naturally expects g νν g bb , g µµ in the absence of large mixings or large hierarchies of U(1) charges. These parameters are then constrained by flavour physics, multiple low-energy precision measurements, colliders, and dark matter detection experiments. In the following sections we work out these constraints, and identify the regions of the parameter space where both the B-anomalies and the dark matter relic abundance can be explained without conflicting any existing experimental data. We note that Z models simultaneously addressing the B-anomalies and dark matter have been previously discussed in refs. [46][47][48][49][50][51][52]. In particular, ref. [49] performed a detailed analysis of collider, precision, dark matter constraints in a similar model based on gauged L µ − L τ symmetry. The main practical difference between our setup and that model is the presence of Z couplings to b-quarks in eq. (2.28), which affects the LHC phenomenology as well as direct and indirect detection signals.

R K ( * ) anomalies and flavour constraints
In this section we review and update the constraints on the parameter space of Z models motivated by the current B-meson anomalies. One possible explanation of the R K and R K * measurements in LHCb is that the low-energy Lagrangian below the weak scale contains an additional contribution to the effective 4-fermion operator with left-handed muon, b-quark, and s-quark fields: Above, the numerical value of the effective coefficient corresponds to the best fit quoted in ref. [7]. In our model, this operator arises from tree-level Z exchange and the analogous operator with µ L replaced by e L does not appear due to vanishing charged lepton mixing. We can express the coefficient G bsµ as function of the couplings in eq. (2.28), Together, eqs. (3.1) and (3.2) imply the constraint on the parameters g bb , g µµ and M Z : There are additional constraints on these parameters coming from flavour physics and low-energy precision measurements. In the following we determine the region of the parameter space where the R K ( * ) anomalies can be explained without conflicting other experimental data.
B s − B s mixing. The Z coupling to bs leads to an additional tree-level contribution to B s − B s mixing. Low-energy observables are affected by the effective operator arising from integrating out the Z at tree level:

JHEP08(2018)061
Such a new contribution is highly constrained by the measurements of the mass difference ∆M s of neutral B s mesons. In this paper we follow the recent analysis of ref. [56] which, using updated lattice results, obtains a stronger bound on G bs : The resulting constraints in the (g µµ , g bb ) plane are shown as the light blue region in figure 1. The updated constraint is particularly strong for the models that generate a strictly positive G bs [56] (as is the case in Z models) due to the ∼ 1.8σ discrepancy between the measured ∆M s and the updated SM predictions which favors G bs < 0. As a consequence, Z models explaining the B-meson anomalies required M Z 1 TeV, assuming weak coupling g µµ 1. For easy reference, we also show the B s mixing constraints based on the previous SM determination of ∆M s [57], − Neutrino trident. The Z coupling to left-handed muons leads to a new tree-level contribution to the effective 4-lepton interaction This operator is constrained by the trident production ν µ γ * → ν µ µ + µ − [58][59][60]. Using the results of the global fit in ref. [61], the bound on the effective coefficient is given by The limits in the (g µµ , g bb ) plane are shown as the orange region in figure 1. Since the trident constraints probe much lower scales than the B s mixing, a much larger Z coupling to muons is allowed, g µµ 1 for a heavy enough Z . Nevertheless, together with the B s mixing constraints, the trident leaves only a narrow sliver of the parameter space that could address the B meson anomalies.
LHC searches. Further constraints on our model come from collider searches. For light Z masses, the LHC measurements of the Z decays to four muons, with the second muon pair produced in the SM via a virtual photon [62,63], pp → Z → 4µ, sets relevant constraints in the low mass region of Z models, 5 M Z 70 GeV. The Z → 4µ constraints on the magnitude of the Z coupling to muons were analyzed in refs. [17,49,60]. Projecting these results onto our model, the excluded parameter space is marked as the pink regions in figure 1 and in the upper-left panel of figure 3. All in all, the Z → 4µ constraint is non-trivial but for any Z mass it always leaves some available parameter space to explain the B-meson anomalies.
For a heavier Z , the strongest constraints comes from LHC dimuon resonance searches, pp → Z → µ + µ − , see also [39]. In our model the Z is dominantly produced at the LHC through its couplings to bottom quarks, bb → Z . The cross section σ(pp → Z ) from bb collisions is taken from figure 3 of ref. [64]. The contribution of bottom-strange collisions, JHEP08(2018)061 Figure 1. The parameter space in the (g µµ , g bb ) plane compatible with R K ( * ) anomalies and flavour constraints (white). The Z mass varies over the plane, with a unique Z mass for each point in the plane as determined by eq. (3.3). We show the recent B s mixing constraints (light blue), and the trident bounds (orange); for reference we also display the previous weaker B s mixing bounds (dark blue). The green, red, purple and black lines correspond to M Z = 10, 100, 1000, 10000 GeV respectively.
which is subleading in our model, is estimated using Madgraph [65]. The Z boson can subsequently decay into muons, muon neutrinos, bottom or strange quarks, and also into top quarks and dark matter when kinematically allowed. The partial decay widths are given by from which we calculate Br(Z → µµ) analytically. Then σ(pp → Z → µµ) is estimated using the narrow-width approximation, and compared with the limits from the recent dimuon resonance search by ATLAS [66], which allows us to constrain Z mases between 150 GeV and 5 TeV. We verified that the analogous Tevatron analyses give weaker constraints, also in the low mass regime. Figure 3 shows the ATLAS constraints for specific Z masses (200, 500 and 1000 GeV) with dark matter couplings set to zero and arbitrary (g µµ , g bb ) couplings. Figure 1 shows the same limits for the Z mass fixed in function of (g µµ , g bb ) by

JHEP08(2018)061
the condition in eq. (3.3). We conclude that in the parameter space of our model relevant for explaining the B-meson anomalies the ATLAS dimuon limits are always weaker that the new B s mixing constraints.
Constraints from lepton-flavour violation. So far we were assuming zero mixing in the charged-lepton sector. It is interesting to discuss the constraints resulting from relaxing that assumption. In particular, for a non-vanishing mixing angle between charged leptons of the second and first generations (V eL ) 21 = 0, a non-diagonal Z coupling to left-handed muons and electrons would be present which could generate an additional contribution to the transition µ → eγ whose partial decay width can be estimated, according to [40], as where F (x) is a loop function, as defined in [40], whose limit for m Z m µ is lim x→0 F (x) = 2/3. The branching ratio of µ → eγ is severely constrained by the MEG experiment [67] which set the bound BR(µ → eγ) ≤ 4.2 × 10 −13 at 90%CL. An analytical approximation of this branching ratio is given by implying that µ → eγ is expected to set a stronger constraint than the neutrino trident production for values of the mixing angle |V eL | 21 10 −4 as represented in figure 2, while |V eL | 21 10 −3 would rule out the entire parameter space. As a result, in the viable parameter space of our setup, the mixing angle |V eL | 21 is expected to be |V eL | 21 10 −4 . Similarly, the experimental limit on the lepton-flavour-violating of the tau lepton into 3 muons, Br(τ → 3µ) ≤ 2 × 10 −8 [68], constrains the mixing angle between charged leptons of the second and third generation (V eL ) 32 : . This is stronger than the trident bound in eq. (3.7) for (V eL ) 32 3 × 10 −4 , while (V eL ) 32 3 × 10 −3 would rule out the entire parameter space.
Other constraints. Finally we comment on other precision observables which yield subleading constraints on our model. The contribution of Z to the muon magnetic moment is given by (3.12) The measured discrepancy of the muon magnetic moment is ∆ µ g−2 = (290±90)×10 −11 [69]. This sets weaker limits on the ratio g µµ /M Z than the trident production.

Dark matter
Our model comprises a fourth neutrino (ν 4 ) which possesses all the properties of a viable dark matter candidate. Indeed, ν 4 is an electrically neutral particle interacting weakly with the SM sector through an exchange of Z . Furthermore, our charge assignments under the local symmetries forbid any mixing with other fields such as the SM neutrinos. Therefore, conservation of fermion number in the dark sector can be effectively seen as a Z 2 symmetry forbidding the dark matter from decaying and as a consequence ensuring its stability. In this section we discuss in some detail the generation of the dark matter relic density and show the constraints from indirect and direct dark matter searches. For brevity, in the following the dark matter candidate is simply denoted as ν.  Figure 3. Bounds on the parameter space in the (g µµ , g bb ) plane for fixed Z masses: 50, 200, 500 and 1000 GeV, as indicated on each panel. The red bands explain R K ( * ) at 1σ. The blue and orange areas show the B s − B s mixing [56] and neutrino trident [61] 2σ exclusions, respectively. For low Z masses we have additional constraints from Z → 4µ as shown in pink. The ATLAS limits [66] from dimuon resonance searches for 36 fb −1 luminosity are given in purple for larger Z masses.

Relic abundance
Our dark matter candidate is a weakly interacting massive particle (WIMP) whose relic abundance can be generated via the well studied freeze-out scenario [72,73]. We will fit the parameters to reproduce the present dark matter density measured by the Planck collaboration: Ω DM h 2 = 0.1198 ± 0.0015 [74]. In the WIMP scenario, Ω DM is determined by the thermally averaged annihilation cross section σv which can be expressed as [75] σv = 1 where x = m ν /T , K 1,2 (x) are the modified Bessel functions of the second kind, and σ s is the dark matter annihilation cross section at the centre-of-mass energy squared s. When σv is approximately independent of the temperature the relic density is related to it as

JHEP08(2018)061
In our model the dark matter particles can annihilate to 3μ µ,ν µ ν µ ,bb, and possibly tō tt, Z Z , if kinematically accessible: One can derive an analytical approximation of σv by expanding it in powers of x −1 around the typical freeze-out temperature x F ∼ 23. Away from the pole and thresholds, each component of σv can be approximated by the s-wave expression: where c ψ is a color factor. One can see that the annihilation cross section grows as m 2 ν for small dark matter masses, and evolves as m −2 ν for large dark matter masses. Therefore, for fixed couplings and M Z , there are typically two possible values of m ν reproducing σv thermal , as illustrated in figure 4. For small couplings, g 0.1, the annihilation cross section is substantially lower than the thermal one except in the pole region, and the two solutions approach m ν ∼ M Z /2. As demonstrated in [76], the presence of a pole in the annihilation cross section may invalidate the 1/x expansion. In such a case one cannot use eqs. (4.4) and instead one has to rely on numerical evaluations using eq. (4.1). In order to explore the complete available parameter space, we compute the relic density and σv numerically using the package micrOMEGAs [77] after implementing the model in FeynRules [78]. For higher values of the couplings, g 1, the correct relic density can be achieved away from the pole region where eqs. (4.4) are adequate.

Indirect detection constraints
In the WIMP framework, dark matter annihilations to SM states occurring inside large astrophysical structures such as the galactic center, dwarf spheroidal (dSphs) galaxies or galaxy clusters might be relatively frequent at the present time. This could lead to indirect dark matter observation by detecting the by-products of these annihilations in high-energy cosmic rays [83]. The absence of any significant signal so far allowed experimental collaborations to derive upper limits on σv as a function of the dark matter mass, assuming a particular dark matter density distribution for these astrophysical structures. In our model σv is approximately velocity independent in the non-relativistic limit, therefore the same value of order σv thermal required to match the relic density is also relevant for indirect detection. The two annihilation channels most relevant for indirect detection are νν → bb and νν → µ + µ − . In the parameter space where νν → bb (and possibly to tt) dominates, the best current limits on σv are derived by the Fermi-LAT collaboration from a combined analysis of 15 Milky Way dSphs and excludes dark matter masses JHEP08(2018)061 m DM 100 GeV [80], assuming the Navarro-Frenk-White profile [84]. For larger dark matter masses stronger constraints on the same annihilation channel come from the HESS experiment [79], however the typical limits are σv 10 −25 cm 3 s −1 and therefore cross sections of order σv thermal are not probed. In the future, sensitivity of the Cherenkov Telescope Array (CTA) might be sufficient to probe annihilation the thermal cross section for m DM 100 GeV [81,[85][86][87][88]. The current and future constraints in the bb annihilation channel are illustrated in figure 4, where we also show predictions of our model for two particular points in the parameter space.
As can be seen in figure 1, given the newer B s mixing constraints the allowed parameter space has g µµ g bb , and therefore annihilation into µ + µ − (and the corresponding neutrinos) dominates. In such a case, the indirect detection limits on σv are substantially weaker, such that the thermal annihilation cross section is allowed for dark matter masses above a few GeV [80]. For this reason, the indirect limits are not relevant in most of the interesting parameter space of our model. However, for small dark matter masses m ν ∼ GeV annihilation into leptons at redshift z ∼ 1000 can be constrained by CMB spectrum observations, as it could modify the ionization history. For the thermal annihilation cross section, the Planck collaboration constraints on CMB spectrum distortions exclude dark matter masses below m ν 10 GeV [82], as illustrated in figure 4. We note that annihilation into leptonic final states can be relevant for experiments such as AMS-02 measuring cosmic-ray positrons and electrons, from which several studies have obtained JHEP08(2018)061 strong constraints on σv [89][90][91][92]. However these constraints are subject to strong dependence on the propagation model and uncertainties regarding cosmic-ray propagation in the interstellar medium, and for this reason we do not include them in the following. All in all, in our model dark matter masses 10 GeV are excluded by the Planck collaboration results. Moreover, in the parameter space where annihilation into bb dominates, dark matter masses below 100 GeV are excluded by the Fermi-LAT results, although that parameter space is also disfavored by the recent B s mixing constraints.

Direct searches
Direct detection (DD) of dark matter has proved to be extremely useful to constrain WIMP scenarios. These experiments aim at observing the recoil energy due to dark matter particles present in the Milky Way halo scattering off nuclei of a detector material. Their sensitivity has improved by several orders of magnitude during the past decade, and currently the xenon-based experiments LUX [93], PandaX [94] and XENON1T [95] probe dark matter spin-independent (SI) scattering cross section of the order of σ SI 10 −45 cm 2 for dark matter masses of the order of 100 GeV.
In our set-up, integrating out the tree-level Z exchange between dark matter and the SM leads to the following effective operators at the scale µ M Z : where again one can neglect the effective coupling to bs. Below the scale M Z , vector-like dark matter couplings to light quarks are induced via renormalization group (RG) running: 1,f (µ)qγ α qνγ α ν . The complete RG equations can be found e.g. in [96]; schematically, one has C Other tensor structures beyond that in eq. (4.6) also appear but they give subleading effects in direct detection. Finally, at µ 2 GeV the couplings in eq. (4.6) can be mapped to momentum-and spin-independent non-relativistic interactions of dark matter with protons and neutrons: where c p 1 = 2C 1,u + C 1,dχ | µ 2GeV and c n 1 = C 1,u + 2C 1,d | µ 2GeV . We evaluate numerically the one-loop RG evolution of effective couplings. To this end, above m Z we use the RunDM package [96][97][98], while running below m Z and the coefficients c  Figure 5. Direct and indirect detection constraints on the parameter space: the orange line represents the model points featuring the correct dark matter relic density and the appropriate Wilson coefficient explaining the R K discrepancy. The gray region shows the parameter space excluded by the XENON1T experiment [95] and the green region represents the parameter space not excluded by direct detection experiments but in tension with the Planck collaboration [82] results.
while for M Z = 1 TeV: The coupling to neutrons vanishes within our approximations when M Z ≤ m Z . For M Z > m Z a non-zero c n 1 can be generated, and is dominated by the top Yukawa contributions to the RG running. The dark matter-nucleon spin-independent cross section can be straightforwardly derived from L eff,NR : (4.10) To compare with experimental bounds, which typically assume equal cross section on protons and neutrons, for a target nucleus with Z protons and A − Z neutrons we introduce the averaged cross section In the allowed parameter space relevant for the B-meson anomalies we have g bb g µµ . Assuming that hierarchy, and also m p m χ , for xenon targets an approximate expression for the averaged cross section reads (4.12) In figure 5 we depicted the values of the g νν coupling satisfying the requirement of having the observed dark matter density as well at the correct value of the couplings g µµ and g bb explaining the R K discrepancy. The left panel of that figure illustrates that, for low M Z , the XENON1T collaboration excludes dark matter masses away from the Z pole but still JHEP08(2018)061  allows for low dark matter masses m ν 10 GeV. However, as discussed in the previous subsection, such low masses are excluded by the indirect Planck constraints, therefore the complementarity of direct and indirect detection searches indicates that the dark matter mass has to be close to the pole m ν ∼ M Z /2. For larger M Z , dark matter masses away from the pole region are allowed, see the right panel of figure 5.

Discussion and conclusion
Our main results are shown in figures 6-10 which show for which parameters our model can address the B-meson anomalies while satisfying all experimental and cosmological constraints. As discussed below eq. (2.28), the relevant parameter space is effectively fivedimensional, and spanned by the Z couplings to dark matter (g νν ), muons (g µµ ), and b quarks (g bb ), and by the masses of dark matter (m ν ) and the Z vector messenger (M Z ). We display it in the {g µµ ,g νν } plane for several representative values of M Z . For each g µµ and M Z , g bb is fixed according to eq. (3.3) to the best fit value reproducing the R K ( * ) measurements. Then m ν is fixed by the requirement of reproducing the correct relic abundance of dark matter. There are typically two distinct solutions for m ν satisfying σv = σv thermal , therefore for each M Z in the left (right) panel we display the solutions with m ν > M Z /2 (m ν < M Z /2). These solutions are color coded in figures 6-10, from JHEP08(2018)061   smaller (blue) to larger (red) m ν . The white regions are where we find no parameters choice to tune the annihilation cross section to the thermal value. The continuously gray-shaded regions are excluded by direct detection, indirect detection, B s mixing, dimuon searches at the LHC, Z decay to four muons and/or muon trident constraints. However, we choose not to shade the region excluded by the recent update of the B s mixing constraints in [56], and instead represent those by a dashed blue line labeled "B s mixing 2017". The region JHEP08(2018)061 represented on the left of this line is excluded by these constraints. For any value of the Z mass in the considered range there exists a range of parameters reproducing the R K ( * ) anomalies and the relic abundance, and passing all experimental constraints to date. However, for lower M Z the allowed region corresponds to g νν g µµ , once the direct (XENON1T) and indirect (Planck) detection constraints together with updated B s mixing constraints are taken into account. In our model the Z coupling to muons is suppressed by a mixing angle between the SM 2nd generation lepton doublet and the 4th generation vectorlike lepton doublet, and thus we expect g νν g µµ . Conversely, g νν g µµ is unnatural and would require a large hierarchy between the corresponding U(1) charges, q ν 4 q L 4 . On the other hand, for 300 GeV M Z 1 TeV we find some allowed parameter space where g νν is a factor of few larger than g µµ , which is plausible. Further increasing M Z requires a sizable Z coupling to muons in order to address the B-meson anomalies, g µµ 1.
Then we are forced back into the unnatural g νν ∼ g µµ region, simply due to perturbativity constraints on g νν rather than some experimental bounds.
To summarize, assuming our model is indeed the correct explanation of the observed R K ( * ) anomalies and dark matter relic abundance, our analysis hints at a particular corner of the parameter space where 300 GeV M Z 1 TeV, m ν 1 TeV, g νν 1, g bb ∼ 0.1g µµ and 0.1 g µµ 1. This viable space implies large mixing with the vector-like fermions to avoid the gauge coupling g getting into the non-perturbative limit, since g µµ = g q L 4 (s L 24 ) 2 . The mixing angles are proportional to the VEVs of the scalar fields, φ ψ , while inversely proportional to the mass of the vector-like fermions. Furthermore, the mass of the Z is generated by the VEVs of the scalar fields, so that M Z ∼ g φ ψ , which sets the scale of the U(1) breaking not far from the TeV scale, φ ψ ∼ TeV. This set an upper limit in the vector-like fermions at around this scale to get the necessary large mixing. This limit is far from the current heavy charged mass bounds which sets M L 4 100 GeV [71,100,101]. Incidentally, that parameter space can be probed by several distinct methods. First of all, the allowed window can be further squeezed by better precision measurements of the trident ν µ N → µ + µ − ν µ N process, and by improving the theoretical precision of the SM prediction for the B s meson mass difference. The above statement is in fact valid for all models where the B-anomalies are addressed by a tree-level Z exchange. What is more specific to models where the Z interactions with the SM fermions originates from mixing of the latter with vector-like fermions is a non-vanishing Z coupling not only to muons but also to b-quarks. This results in a non-negligible rate of the partonic process bb → Z → µ + µ − which can be probed by dimuon resonance searches at the LHC. In fact, the preferred M Z range is where the LHC sensitivity is optimal. Targeted searches for b-quark-collision initiated process (rather than recast of generic dimuon searches) could lead to a discovery signal in the near future, or to better constraints that are more stringent than the B s mixing one. Finally, the preferred range of dark matter masses and couplings can be probed by direct detection experiments, such that the improvements of one or two orders of magnitude in sensitivity in the next years, which is expected to be achieved by the LZ [102], DARWIN [103] and DarkSide-20k [104] experiments. As illustrated in figure 11, these future improvements should exclude the remaining most natural parameter space of our model. JHEP08(2018)061 Figure 11. Projection of future constraints on the parameter space of our model for M Z = 500 GeV. The current ATLAS dimuon limits [66] are scaled with integrated luminosity to L = 200 fb −1 (ATLAS RUN-2) and L = 3000 fb −1 (ATLAS HIGH L). Future direct detection limits (FUTURE DD) assume that the current XENON1T [95] constraints on the DM-nucleon scattering cross section are improved by two orders of magnitude, which roughly corresponds to the projected sensitivity of LZ [102], DARWIN [103], and Darkside-20k [104] experiments.