Scalar dark matter behind $b \to s \mu \mu$ anomaly

We construct a scalar dark matter model with $U(1)_{L_\mu-L_\tau}$ symmetry in which the dark matter interacts with the quark flavours, allowing lepton non-universal $b \to s \ell \bar{\ell}$ decays. The model can solve $b \to s \mu \mu$ ($R_{K^{(*)}}$) anomaly and accommodate the relic abundance of dark matter simultaneously while satisfying the constraints from other low energy flavour experiments and direct detection experiments of dark matter. The new fields include vector-like heavy quarks $U$ and $D$, $U(1)_{L_\mu-L_\tau}$ breaking scalar $S$, as well as the dark matter candidate $X_I$ and its heavy partner $X_R$. To explain both $b \to s \mu \mu$ anomaly and the dark matter, {\it i)} large mass difference between $X_R$ and $X_I$ is required, {\it ii)} electroweak scale dark matter and heavy quarks are favoured, {\it iii)} not only electroweak scale but ${\cal O}(10)$ TeV dark gauge boson $Z'$ and $X_R$ are allowed.


Introduction
The flavour changing neutral current (FCNC) processes are known to be sensitive to new physics (NP) because they first occur at loop level in the standard model (SM) and therefore are sensitive to heavy physics in the loop. The NP scale they can probe is usually much higher than the scale the LHC can produce. And these indirect searches for NP are complementary to the collider searches. Among many FCNC processes, the b → s transition has been drawing much interest for the last several years because of anomalies in B → K ( * ) µµ and B s → φµµ decays.
In particular SM predictions on the ratio of branching fractions are close to unity, signifying the lepton flavor universality (LFU) in the SM. However, the measurements at the LHCb for K [1] and K * [2] are lower than unity at 2.3 − 2.6σ level. Because the ratio (1.1) is free from hadronic uncertainty, it would be a clear sign for NP, if this violation of LFU persists in future experiments. Including other observables, such as an angular observable in B → K * µ + µ − and branching fraction of B s → φµ + µ − and Λ b → Λµ + µ − , the deviations from the SM predictions increase as large as about 5σ [3][4][5][6], which we will call b → sµµ anomaly. At m b scale the b → s transition is described by the effective weak Hamiltonian  [7,8]. The results from global fitting analyses [3][4][5][6] show that sizable NP contributions to C µ 9 and/or C µ 10 can explain the b → sµµ anomaly.
In this paper we consider a NP model with C µ,NP 10 = 0, in which case the best fit value for C µ,NP 9 is [4] C µ,NP 9 = −1.11 ± 0.17, (1.4) with a SM pull of 5.8σ. In addition to the SM gauge groups we introduce a new gauge symmetry U(1) Lµ−Lτ under which the 2nd (3rd) generation leptons are charged with +1(−1). It is known that the theory is anomaly-free even without extending the SM fermion contents. Since the U(1) Lµ−Lτ gauge boson Z couples to muon, it can make a contribution to the anomalous magnetic moment of muon (g − 2) µ [9,10]. But the Z should be very light ( 400 MeV) to fully accommodate the discrepancy between the experiments and the SM predictions in the (g − 2) µ [11]. The model can also be extended to accommodate neutrino data [12][13][14][15]. In ref. [16] we introduced a fermion dark matter (DM) model whose stability is originated from U(1) Lµ−Lτ symmetry [17]. The model can also explain b → sµµ anomaly by introducing SU(2) L -doublet scalar field, and we showed that there is a strong interplay between the DM and B-physics phenomenology [18][19][20][21][22][23][24][25][26][27][28][29][30][31][32][33]. In this paper we consider a "spin-flipped" version of the model in ref. [16]. We introduce two complex scalar fields S and X: S breaks the U(1) Lµ−Lτ symmetry spontaneously by developing vacuum expectation value (VEV) S , while the lighter component X I is stable by the remnant discrete Z 2 symmetry after U(1) Lµ−Lτ symmetry is broken spontaneously and become a DM candidate. To explain the b → sµµ anomaly as well in this model, we introduce a vector-like quark Q which can mediate quark couplings to Z boson. We will study the solution of the b → sµµ anomaly and the DM phenomenology in this model. This paper is organized as follows. In section 2, we introduce the model and calculate the new particle mass spectra. In section 3 we calculate NP contribution to b → sµµ, and consider low energy constraints including C ,NP

9
, ∆m s in B s −B s mixing, B → K ( * ) νν, b → sγ, the anomalous magnetic moment of muon a µ , and the loop-induced effective Zbb coupling. In section 4 we consider dark matter phenomenology. Finally we conclude in section 6. Loop functions are collected in appendix A.

The model
We introduce a scalar dark matter candidate X and a scalar boson S which gives a mass to U(1) Lµ−Lτ gauge boson Z after the symmetry is broken down spontaneously by the

JHEP05(2019)104
New fermion New scalars Q X S VEV of S. To couple the Z gauge boson to the quarks we also introduce a vector-like where D µ is the covariant derivative, i(= 1, · · · , 3) is the quark-generation index, and is the field strength tensor. The scalar potential is in the form, (2. 2) The trilinear µ term allows a remnant discrete Z 2 symmetry after S gets VEV and U(1) Lµ−Lτ is spontaneously broken. This local Z 2 symmetry [34] stabilises the DM candidate, which we assume the lighter component of X. The kinetic mixing angle χ is strongly constrained to a level of O(10 −3 ) by the DM direct search experiments [35]. The nonvanishing χ does not help solving b → sµµ anomaly because the SM gauge bosons allow only LFU couplings. In this paper we neglect this term for simplicity by setting χ ≡ 0. We note that the fermion Q which has the same SM quantum numbers with the left-handed quark doublets is vector-like under both U(1) Lµ−Lτ and the SM gauge groups. We now consider the particle spectra and identify the DM candidate. After H and S get VEVs, v H and v S , the µ term makes the complex scalar X split into two real scalar fields X R,I defined by with masses-squared

JHEP05(2019)104
Assuming µ > 0, the X I which is the lightest Z 2 -odd neutral field is identified as the DM candidate. The other Z 2 odd fields after S gets VEV are X R and Q. The remaining particles including the SM fields are Z 2 -even. We take m R,I , λ HX , λ SX as free parameters, then we can write the parameters m 2 X and µ in the Lagrangian as After S gets VEV, the U(1) Lµ−Lτ gauge boson Z obtains mass, where g Z is the gauge coupling constant of the U(1) Lµ−Lτ group. The vector-like quark Q does not mix with the ordinary quark with the same SM quantum numbers because it is Z 2 -odd while the SM counterparts are Z 2 -even. So it is already in the mass eigenstates with mass M Q at tree level, though the mass-splitting can be generated at loop-level. This also distinguishes our model from the models in [19,36].
In the unitary gauge we decompose the SM Higgs H and the dark scalar S as The stationary condition at the vacuum gives conditions Using the above conditions it is straightforward to obtain the scalar mass-squared matrix  where c H (s H ) is an abbreviation of cos α H (sin α H ). We require λ HS > −2 √ λ H λ S to stabilise the scalar potential at the electroweak (EW) scale.
In (2.1) we assume that the down-type quarks are already in the mass basis and that the flavor mixing due to Cabibbo-Kobayashi-Maskawa (CKM) matrix V appears in the upquark sector, i.e. d iL = d iL , u iL = j V * ji u jL where primes represent the mass eigenstates. Then the Yukawa interactions of the type q − Q − X can be written in the form where λ u i ≡ j V ij λ j and λ d i ≡ λ i . We will also use the notation λ d,s,b for λ d i and λ u,c,t for λ u i (i = 1, 2, 3). We will simply set λ 1 = 0 to remove the constraints related to the first generation quarks. Even in this case we see is induced. The induced λ u can generate NP contribution to D 0 −D 0 mixing. However, due to Cabibbo-suppressed contribution to D 0 −D 0 at least by O(λ 2 ) where λ(≈ 0.23) is the Cabibbo angle, the constraint from D 0 −D 0 can be always satisfied once the constraint from B s −B s is imposed [24]. And we do not consider this constraint further. The DM interacts with the SM fields through the Higgs-portal Lagrangian (2.14) In this paper we will set α H = 0 and λ HX = 0 to suppress the stringent constraint from the dark matter direct detection experiments via this Higgs portal interaction [37,38].
Assuming m I , m R , and M D are at the EW scale (≡ M EW ), we can neglect terms proportional to external quark mass squareds, m 2 s(b) /M 2 EW ( 1). In this approximation, it is straightforward to get the effective vertex for diagrams (a) and (b) in figure 1: where we take the dimension of space-time integration d to be d ≡ 4 − 2 for positive infinitesimal . The C's are abbreviations for one-loop three-point functions defined in [39], where k = 0, 12, 00 and i, j = I, R(i = j). We will set m s = m b ≡ 0 in the calculation of the C-functions to be consistent with our approximation m 2 1. The C 00 -functions are divergent while C 0 -and C 12 -functions are finite. The divergence in the C 00 -functions can be isolated as where C 00 | finite is the remaining finite part. Using the relation between the U(1) Lµ−Lτ charges, q Q = −q X , we can show that the sum of the two one-loop effective vertices is finite and given by Now we can attach the external muon line in figure 1 to the Z to get C µ 9 . The full amplitude for b → sµµ transition in figure 1 is given by where p 3 (p 4 ) is outgoing (incoming) muon four-momentum. The term proportional to q µ q ν vanishes becauseū(p 3 ) q u(p 4 ) = 0. Since q ∼ O(m b ) at most, we can set q 2 ≡ 0 in the denominator of Z -propagator. In this case the effective vertex can be written in a simple analytic form: where x I(R) = m 2 I(R) /M 2 D and the loop function k is defined in appendix A. We note V sb → 0, when x R → x I . This can be understood as follows: in the limit m R → m I , the two real scalars X I and X R merge into the original complex scalar X as can be seen from (2.3). In this limit, the subset of the full Lagrangian which contributes the effective vertex iV µ eff (q 2 ), is invariant under local U(1) Lµ−Lτ symmetry. Note that the Z mass term which breaks the U(1) Lµ−Lτ is not included in ∆L. Then the Ward-Takahashi identity dictates V eff (q 2 = 0) = 0, 1 due to the absence of the tree-level Z -exchanging FCNC, leading to V eff (q 2 ) ∝ q 2 to all orders of perturbation theory [40]. Since V eff (q 2 = 0) = 0, we obtain V sb = 0 in the limits q 2 → 0 and m R → m I . Now it is straightforward to get where the prime on the k functions denotes a derivative with respect to the argument and we fixed the U(1) Lµ−Lτ charge of µ − to be +1. A sizable mass splitting between m R and m I is favoured to generate C µ,NP 9 which can explain the b → sµµ anomaly. As a benchmark point in the parameter space, we choose q X = 2, α Z = 0.1, m Z = 700 GeV, λ s λ * b = 0.2, M D = 1 TeV, m I = 900 GeV, and m R = 3 TeV, for which we get which is close to the best fit value in (1.4) to solve the b → sµµ anomaly.

JHEP05(2019)104
m Z = 700 GeV, M D = 1 TeV, q X = 2, α Z = 0.1, and λ s λ * b = 2. As can be seen from (3.9), |C µ,NP 9 | becomes larger as the mass splitting m R − m I increases. Figure 2(b) shows a plot in the same plane but by varying M D = 1, 2, 3 TeV (from below) with M R = 3 TeV and the other parameters the same as figure 2(a). We can see that the effect of the vectorlike D-quark decouples as M D increases. In both cases, smaller m I is favored to obtain larger |C µ,NP 9 |. Photon-or Z-penguin diagrams similar to Z -penguin diagrams in figure 1 but with Z replaced by photon or Z-boson can contribute to C µ,NP

9
. Their couplings to leptons are flavour-universal and they also contribute to C e,NP 9 and C τ,NP 9 with the same value. So we use them as a constraint on the model. The one-loop effective vertices they generate are proportional to q 2 by the same logic used to show V eff (q 2 ) ∝ q 2 above. Here the conserved U(1) symmetries are the U(1)-electromagnetism, U(1) em , for photon vertex, and the neutral current part of SU(2) L , U(1) Z , for Z-boson vertex. Since these symmetries are conserved whether U(1) Lµ−Lτ is conserved or not, the argument applies even when m I = m I . If we attach the external muon lines, the q 2 in the photon-vertex cancels q 2 in the photon propagator, whereas the one in the Z-vertex does not. As a consequence, the Z-penguin contribution is negligible because it is proportional to We obtain the photon penguin contribution to be GeV, and m R = 3 TeV, we get which is about three orders of magnitude smaller than the Z contribution to C µ,NP 9 in (3.10). And we can neglect the photon-and Z-penguin contributions. Now we consider other constraints on the model parameters. It turns out that the value |λ s λ * b | is the most strongly constrained by the measurements of the mass difference ∆m s for B s −B s mixing. Figure 3 shows one-loop box diagrams for B s −B s mixing. The arrows represent color flow. The lower two diagrams with crossed scalar lines exist because X I and X R are real scalars. Our model where new particles couple only to the left-handed quarks contributes to the same effective operator with the one in the SM, (3.13) The Wilson coefficient C 1 can be decomposed into the SM and the NP contributions (3.14) The SM contribution at the electroweak scale is obtained by box diagrams with W -boson and t-quark running inside the loop:  where x t = m 2 t /m 2 W ≈ 4.64 and the loop function S 0 (x t ) can be found in [41]. The NP diagrams shown in figure 3 give We note that this result is non-vanishing, different from a single real DM contribution which vanishes [24]. The non-zero term arises from the diagrams with X R and X I at the same time. The measurement of the mass difference in the B s −B s system gives a constraint on the value of C NP 1 : −2.1 × 10 −11 ≤ C NP 1 ≤ 0.6 × 10 −11 (GeV −2 ), (3.17) at 2σ confidence level [24]. For the benchmark point m I = 900 GeV, M D = 1 TeV and m R = 3 TeV, we get which is about an order of magnitude smaller than the SM prediction. This point satisfies the constraint (3.17).  For this process the SM prediction has been calculated up to NNLO QCD corrections [43], which predict, The NP contribution to B(B → X s γ) can be obtained by calculating the Wilson coefficients C 7γ,8g from the diagrams in figure 5: where e D = −1/3 is the electric charge of the vector-like down-type quark D and x i = m 2 i /M 2 D (i = I, R). The loop-fuction J 1 (x) is listed in the appendix A. From the prediction including NP contribution to C 7γ(8g) [43], (3.19) and (3.20), we obtain the constraint which is about two orders of magnitude less than the current bound (3.22). Figure 6 shows a contour plot for the combination C 7γ + 0.24 C 8g with contour lines −10 −5 , −5 × 10 −5 , and −10 −4 in the (m I , M D )-plane for λ s λ * b = 0.2 and three different values of m R . The green solid (yellow dashed, magenta dot-dashed) lines correspond to m R = 1(2, 3) TeV. We can see C 7γ + 0.24 C 8g is less sensitive to m R than C NP 1 of B s −B s is. The entire region considered is allowed by B exp (B → X s γ).
The NP diagrams for semi-leptonic decay B → K ( * ) νν is obtained when the external muon lines are replaced with neutrino lines in figure 1. The effective Hamiltonian for the decay is where (3.25)

JHEP05(2019)104
We obtain the non-vanishing coefficients We note that the diagram with Z replaced by Z vanishes in the q 2 → 0 limit, showing that the Z contribution is dominant. The current experimental bounds on the ratios of branching fractions In our model these ratios are predicted to be showing the deviation from the SM is very small partly due to the cancellation of the interference terms. The gauged U(1) Lµ−Lτ model is well known to generate the sizable anomalous magnetic moment of muon, a µ = (g − 2) µ /2 via the Z -exchanging one-loop diagram [9]. The Z contribution can explain the long-standing discrepancy between the experimental measurements [44] and the SM predictions [45]: The effective Hamiltonian for a µ is The NP contribution at one-loop level is calculated to be  Figure 7. The constraint from neutrino trident production and muon g − 2 in (m Z , α Z ) plane. The grey region is disfavoured by the neutrino trident production experiments at 2σ. The region between the two green lines is favoured by the current discrepancy ∆a µ at 2σ, but it is excluded by the neutrino trident production experiments. The red "X" mark represents the benchmark point m Z = 700 GeV and α Z = 0.1.
For the benchmark point m Z = 700 GeV, α Z = 0.1, we get which is consistent with (3.31) within 3σ. In the minimal U(1) Lµ−Lτ model the region in the (m Z , α Z ) plane which can explain ∆a µ at 2σ level is excluded by the bound from the measurement of neutrino trident production, ν µ N → ν µ N µ + µ − , when m Z 400 MeV [11]. Figure 7 shows the constraint from neutrino trident production and muon g−2 in (m Z , α Z ) plane. The grey region is disfavoured by the neutrino trident production experiments at 2σ. The region between the two green lines is favoured by the current discrepancy ∆a µ at 2σ, but it is excluded by the neutrino trident production experiments. The red "X" mark represents the benchmark point m Z = 700 GeV and α Z = 0.1.
The new particles in the model also generates one-loop effective Zff -vertex (f = s, b). Since Zbb vertex has been more precisely determined by the LEP experiment, we consider the constraint only from Zbb. The Zbb vertex is written in the form, with a correlation coefficient of +0.9. In our model the NP contributions to g b L is obtained to be, 3 where x I(R) = m 2 I(R) /M 2 D and we can take q 2 = m 2 Z . The loop function Q 1 (x) is listed in (A.5). We notice that the loop function is the same with the one for the photon penguin diagram of b → s in (3.11). In both cases the gauge bosons couple only to D and the amplitudes are proportional to q 2 by Ward-Takahashi identity as we mentioned above (3.11). So they should be proportional to each other. The above result can be compared with (3.37). In figure 8 we show error ellipses at 1, 2, and 3σ confidence level in (δg b R , δg b L ) plane. The vertical red line segment is obtained by randomly scanning g b,NP L in the ranges 10 −3 < λ b < 1, 10 GeV < m I < 3 TeV, (3.39) The model predicts δg b,NP L in the range (−4.9 × 10 −4 , 0) and satisfies (3.37) at 3σ level. The new particle searches at the LHC can also constrain the model. For example, new coloured-scalars D or U can be pair-produced via pp → DD(UŪ ) at the LHC if their masses are within the LHC reach. These production processes are similar to those considered in [47,48] where they were analysed in detail. Roughly M D(U ) 1 TeV are excluded. And we impose M D ≥ 1 TeV.

JHEP05(2019)104 4 The dark matter
In this section we identify the main channel and the favoured parameter region to give the observed DM relic density, Ω DM h 2 = 0.1199 ± 0.0022 [49]. We assume the weakly interacting massive particle (WIMP), X I , whose mass is at the electroweak scale, is the candidate for a cold dark matter (CDM) and constitute the whole dark matter component in the universe. In addition we assume the DM relic came from the thermal freeze-out mechanism. In this mechanism, when they are at the initial equilibrium state for the high temperature, T ∼ m I , the DM particles whose number density is similar to that of the photon are overabundant. The DM number density becomes reduced by (co)annihilations until their rates are smaller than the Hubble expansion rate, when it freezes out typically near T ∼ m I /25 [50]. Then the relic density is roughly related with the annihilation cross section at freeze-out temperature as where v is the relative velocity between the DM particles. Before studying DM phenomenology, we can get insight by comparing our model with the minimal "scalar singlet dark matter" model with Z 2 symmetry [51]. The scalar potential in the minimal model has terms The DM mass m D is obtained by m 2 The DM annihilation occurs through DD → h → SM SM or DD → hh. Both processes are controlled by the Higgs portal coupling λ hD , which is strongly restricted by the direct detection experiments [52][53][54]. As a consequence the model is strongly constrained, ruling out m D 1 TeV region as a single-component DM [55].
In our model, however, there are many model parameters involved in the DM-Higgs couplings as can be seen in (2.14), which allows the direct detection constraint on the Higgs portal interaction to be lifted by setting α H = λ HX = 0 to remove H 1 X 2 I term. Even in that case the heavy Higgs H 2 can still mediate the DM interaction without much affecting the DM scattering off nuclei. There are also dark gauge interaction and dark Yukawa couplings available for DM annihilations. In this paper we consider two processes for DM annihilation which can occur in different regions of parameter space: X I X I → Z Z and X I X I → qq. Barring the Higgs portal X I interaction with H 1 , they are dominant processes. In figure 9 we show representative diagrams for the two annihilation channels. We implemented our model to the micrOMEGAs package [56] to evaluate the DM relic density and direct detection cross section.
We first consider the scenario in which the diagrams of type (a) in figure 9 play a major role. The process X I X I → Z Z dominates the DM annihilations as long as it is kinematically open, α Z is not too small (α Z 10 −6 ), and m H 2 is not much larger than TeV scale (m H 2 5 TeV). In this case we obtain typically σv(X I X I → Z Z ) σv(X I X I → qq ) for the thermal-averaged annihilation cross sections. Given that we set α H = λ HX = JHEP05(2019)104 Figure 9. Representative diagrams for the two annihilation channels: (a) X I X I → Z Z and (b) u, d, s, c, b, t). 0, the process X I X I → Z Z is controlled by the dark Higgs interaction and the dark gauge interaction. The former interaction is given by λ S v S − √ 2µ, and the latter by g Z v S . Both are sensitive to C µ,NP 9 in (3.9). Figure 10 shows contour lines of Ω DM h 2 = 0.1199 in (m I , α Z ) (left panel) and (m I , m R ) (right panel) plane for m Z = 10, 100, 700 GeV. For the other parameters we take the benchmark values: q X = 2, λ s λ * b = 0.2, M D = 1 TeV, and m H 2 = 2 TeV. We set m R = 3 TeV for the left panel and α Z = 0.1 for the right panel. We also fixed α H = 0, λ X = 1, λ HX = 0, and λ SX = 0.2. At this stage we do not impose constraints other than m R > m I . Numerically we have checked that σv(X I X I → Z Z ) is much larger than σv(X I X I → qq ) for the points in figure 10.
The process X I X I → Z Z can occur in the early universe even when m I < m Z if the mass difference is not too large. This can be seen in the steep lines corresponding to For example, for m I = 90 GeV and m Z = 100 GeV, we obtain v ≥ 0.87. The DM should be quite relativistic and the thermally averaged annhilation cross section is Boltzmannsuppressed. When X I X I → Z Z is kinematically open for non-relativistic X I , the process is sensitive to the dark gauge coupling α Z . For our benchmark point it turns out that the H 2 -exchanging s-channel diagram is more important than the X R -exchanging t-channel diagram due to the µ−term in (2.14). This shows that the process is also sensitive to the mass-squared difference, m 2 R − m 2 I , by (2.5). When m I is not close to the resonance region, the s-wave annihilation cross section for the X I X I → H 2 → Z Z channel is in the form When λ SX is not too large, the larger mass squared difference m 2 R − m 2 I and the smaller m 2 Z , i.e. the larger µ, the larger σv is obtained. As m I approaches 1 TeV, it is close to the resonance region m H 2 ≈ 2m I and the cross section increases rapidly, virtually independent of α Z . This explains almost vertical parts of the curves near m I = 1 TeV.
In the right panel of figure 10, we can see that the regions m R ≈ m I also give the correct relic density. This occurs due to the coannihilation processes X R X I → Z H i (qq, ¯ ) [16]. As we saw in (3.9), the NP contribution to C µ,NP is suppressed. And the coannihilation mechanism for the DM relic density is not favoured as a solution to b → sµµ anomaly. This shows a strong interplay between the flavour physics and DM phenomena [16,48,[57][58][59][60][61][62], which will be discussed in the next section in more detail. Once X I X I → Z Z is kinematically open near m I = m Z , it dominates the annihilation processes, which is not so sensitive to m R . Now we consider the parameter space where type (b) diagrams in figure 9 dominate the annihilation cross section. Figure 11 shows the dependence of the dark Yukawa coupling λ b as a function of m I to give constant Ω DM h 2 = 0.1199 for three different values of M D : M D = 1, 2, 3 TeV. For this we take heavy m Z = 10 TeV and m H 2 = 10 TeV to suppress the X I X I → Z Z channel. We set λ s = 0.4. For the other fixed parameters, we take the same values with those for figure 10. The λ b required to give Ω DM h 2 changes sharply near m t /2 and m t . This occurs due to the processes X I X I → ct(ct) and X I X I → tt, respectively. Their s-wave annihilation cross sections are given by where we have neglected the mass of charm quark. Note that both are proportional to m 2 t . The contribution from X I X I → bb, being proportional to m 2 b , is negligible compared to the above two processes. 4 When m I ≈ M D , the near vertical lines are due to the coannihilation processes such as DD(UŪ ) → gg, qq, Zg (q = u, d, s, c, b, t), UD → W + g, and X I D → sg. They are sensitive to the SM SU(3) C gauge coupling and independent of λ b . If we require λ b to be of order 1, to give the correct relic density M D should not be much heavier than O(1) TeV, increasing the prospect of producing D or U at the LHC.

Interplay between the b → sµµ anomaly and the dark matter
Now we investigate whether the parameter space which solves the b → sµµ anomaly can also give the correct relic density for the dark matter. At this stage we impose the low energy flavour constraints discussed in the previous sections. We also consider constraints from the direct detection experiments of dark matter such as LUX [52], PANDA [53], and XENON1T [54]. Figure 12 shows plots for C µ,NP 9 which solves the b → sµµ anomaly at 1σ (dark blue) and 2σ (light blue) in (m I , α Z ) plane. We take m Z = 700, 100, 10 GeV (from the left panel). We fixed q X = 2, λ s λ * b = 0.2, M D = 1 TeV, and m R = 3 TeV. They are superimposed with the constant lines for Ω DM h 2 = 0.1199 shown in figure 10. The grey regions are excluded by the neutrino trident production experiments at 2σ level. We checked that the other low energy experiments do not further constrain the allows regions for the C µ,NP 9 and the relic density. Neither does the direct detection experiments affect the plots in figure 12 because i) the Z couples to the quarks at one-loop level and ii) more importantly only inelastic upward scattering X I q → X R q can occur for Z interaction, which is forbidden kinematically. We can see that the b → sµµ anomaly can be resolved JHEP05(2019)104 for b → sµµ anomaly (blue region) and the relic density of the dark matter (red line) in (m I , m R ) plane. The dark blue, blue, and light blue region represent 1, 2, and 3σ allowed region, respectively. From the left panel we take m Z = 700, 100, 10 GeV. We fixed q X = 2, α Z = 0.1, λ s λ * b = 0.2, and M D = 1 TeV. The other fixed parameters are the same with those for figure 10. The grey region is unphysical because m R < m I . at 1σ for m Z = 700 GeV and the current dark matter can be accommodated at the same time. For smaller Z masses, m Z = 100, 10 GeV, C µ,NP 9 becomes too large and b → sµµ anomaly can be explained only at 2σ level to explain the current relic density. This result shows a strong interplay between low energy B-meson decay experiments and the dark matter physics.
In figure 13 we show the C µ,NP 9 for b → sµµ anomaly (blue region) and the relic density of the dark matter (red line) in (m I , m R ) plane. From the left panel we take m Z = 700, 100, 10 GeV. As in the previous case we fixed q X = 2, α Z = 0.1, λ s λ * b = 0.2, and M D = 1 TeV. The other fixed parameters are the same with those for figure 10. The grey region is unphysical because m R < m I and we exclude it. The regions we considered are not constrained by other observables such as B s −B s mixing, neutrino trident production, or direct detection of dark matter. For each case there is intersection JHEP05(2019)104 region of the required C µ,NP 9 and the correct relic density free from other experimental constraints. The region occurs near the kinematic threshold of X I X I → Z Z . Relatively large m R compared to m I is also required to get sizable C µ,NP

9
. Now we consider the impact of λ couplings on the dark matter and the C µ,NP

9
. For this purpose we suppress the dark gauge contributions to them by decoupling Z as in figure 11. To decouple we assume Z is heavy: for our purpose it is enough to set m Z = 2 TeV as in figure 11. Figure 14 shows the results of this setting in (m I , λ b ) plane. We fixed q X = 2, α Z = 0.1, λ s = 0.4, m R = 20 TeV and M D = 1, 2, 3 TeV (from the left panel). We take the same values with those for figure 11 for the other fixed parameters. The grey region is excluded because the cross section of the DM scattering off the nuclei is too large. For M D = 3 TeV the direct detection constraint disappears completely from the region considered. The region with peach color is excluded by the experimental constraints on ∆m s of B s −B s system. We can see that the direct detection experiments and the B s −B s mixing play a complementary role in excluding the parameter region, although the latter plays more important role in our choice of parameters. Both the relic density and the b → sµµ anomaly can be explained simultaneously for the electroweak scale DM. We notice that the contribution to C µ,NP 9 is not easily decoupled for very heavy m Z and m R due to large mass splitting m R − m I . Eventually as m Z and/or m R becomes even heavier, their impact on C µ,NP 9 will get smaller. To see this decoupling effect we need to resum the large logarithm of log(m R /M D ) and also consider higher loop effects, which is beyond the current analysis.

Conclusions
We considered a new physics model with the U(1) Lµ−Lτ symmetry which has both a dark matter candidate and new flavour changing neutral currents in the quark sector. This opens up a possibility that there may exist a strong interplay between the dark matter and the flavour physics. In particular we showed that we could simultaneously explain the b → sµµ anomaly and the dark matter abundance in our universe. The model has a scalar dark matter candidate X, a SU(2) L -doublet colored fermion Q = (U, D) T , and a dark Higgs S JHEP05(2019)104 with k(1, 1) = 3/2. The loop function J 1 (x) for b → sγ is We have J 1 (1) = 1/24. The loop function for the photon-penguin of b → s ¯ the effective Zbb vertex is We get Q 1 (1) = 1/8.
Note added. After finalizing the manuscript we received a paper considering b → sµµ anomaly in a similar but a different setting [63]. Their model does not have µ term, and as a consequence the dark matter candidate is complex scalar while it is real scalar in our case. Their results show the allowed region is rather restricted compared to ours due to the absence of the µ term.