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

TL;DR

This paper proposes a scalar dark matter model with a specific symmetry that explains the $b o s \mu \mu$ anomaly and matches dark matter relic abundance, while satisfying experimental constraints.

Contribution

It introduces a novel scalar dark matter model with $U(1)_{L_\mu-L_ au}$ symmetry linking dark matter to flavor anomalies and experimental constraints.

Findings

The model can simultaneously explain the $b o s \mu \\mu$ anomaly and dark matter relic abundance.

Large mass differences between dark matter components are necessary.

Electroweak scale and multi-TeV scale particles are compatible with the model.

Abstract

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.

Figures

Tables

Table 1. Table 1: New particles in the model with their quantum numbers under the gauge group S ​ U ​ ( 3 ) C × S ​ U ​ ( 2 ) L × U ​ ( 1 ) Y × U ​ ( 1 ) L μ − L τ 𝑆 𝑈 subscript 3 𝐶 𝑆 𝑈 subscript 2 𝐿 𝑈 subscript 1 𝑌 𝑈 subscript 1 subscript 𝐿 𝜇 subscript 𝐿 𝜏 SU(3)_{C}\times SU(2)_{L}\times U(1)_{Y}\times U(1)_{L_{\mu}-L_{\tau}} .

	New fermion	New scalars
	$Q$	$X$	$S$
$SU{\left(3\right)}_C$	3	1	1
$SU{\left(2\right)}_L$	2	1	1
$U{\left(1\right)}_Y$	$\frac{1}{6}$	$0$	$0$
$U{\left(1\right)}_X$	$q_Q\left(\equiv -q_X\right)$	$q_X$	$q_S\left(\equiv -2q_X\right)$

Equations

$$R_{K^{(\ast )}}=\frac{\mathcal{B}(B\to K^{(\ast )}{\mu }^{+}{\mu }^{-})}{\mathcal{B}(B\to K^{(\ast )}{e}^{+}{e}^{-})},$$

$${\mathcal{H}}_{\mathrm{eff}}$$

$$O_{7\gamma }^{(\prime )}$$

$$O_9^{(\prime )\ell }$$

$$C_9^{\mu ,\mathrm{NP}}$$

$$\mathcal{L}$$

$$-\sum _{i=1}^3({\lambda }_i{\bar{q}}_L^{i}QX+h.c.),$$

$$V$$

$$+{\lambda }_{HS}{H}^{†}H{S}^{†}S+{\lambda }_{HX}{H}^{†}H{X}^{†}X+{\lambda }_{SX}{S}^{†}S{X}^{†}X.$$

$$X\equiv \frac{1}{\sqrt{2}}(X_R+i{X}_I),$$

$$m_{R,I}^2=m_X^2+\frac{1}{2}{\lambda }_{HX}{v}_H^2+\frac{1}{2}{\lambda }_{SX}{v}_S^2\pm \sqrt{2}\mu v_S.$$

$$m_X^2$$

$$\mu$$

$$m_{Z^{\prime }}=g_{Z^{\prime }}|q_S|v_S=2g_{Z^{\prime }}|q_X|v_S,$$

$$\displaystyle H=\left(\begin{array}[]{c}0\\ {1\over\sqrt{2}}(v_{H}+h)\end{array}\right),\quad S={1\over\sqrt{2}}(v_{S}+s).$$

$${\mu }_H^2$$

$${\mu }_S^2$$

$$\left(\begin{array}[]{cc}2\lambda_{H}v_{H}^{2}&\lambda_{HS}v_{H}v_{S}\\ \lambda_{HS}v_{H}v_{S}&2\lambda_{S}v_{S}^{2}\\ \end{array}\right),$$

$$\displaystyle\left(\begin{array}[]{c}h\\ s\end{array}\right)=\left(\begin{array}[]{cc}\cos\alpha_{H}&\sin\alpha_{H}\\ -\sin\alpha_{H}&\cos\alpha_{H}\\ \end{array}\right)\left(\begin{array}[]{c}H_{1}\\ H_{2}\end{array}\right),$$

$${\lambda }_H$$

$${\lambda }_S$$

$${\lambda }_{HS}$$

$$\Delta {\mathcal{L}}_{\mathrm{Yukawa}}$$

$${\lambda }_u=V_{us}{\lambda }_2+V_{ub}{\lambda }_3,$$

$$\mathcal{L}(X_I{X}_I{H}_{1,2})=$$

$$\displaystyle-{1\over 2}\Big{[}\lambda_{HX}v_{H}s_{H}+(\lambda_{SX}v_{S}-\sqrt{2}\mu)c_{H}\Big{]}H_{2}X_{I}^{2}.$$

$$i{V}_{\mu }^{(a)}$$

$$i{V}_{\mu }^{(b)}$$

$$C_k^{(i)}$$

$$C_k^{(ij)}$$

$$C_{00}=\frac{1}{ϵ}-{\gamma }_E+\log 4\pi +C_{00}{|}_{\mathrm{finite}},$$

$$i{V}_{\mathrm{eff}}^{\mu }(q^2)$$

$$=\frac{-i{q}_X{g}_{Z^{\prime }}{\lambda }_s{\lambda }_b^{\ast }}{32{\pi }^2}{\mathcal{V}}_{sb}(q^2,M_D^2,m_I^2,m_R^2)\bar{u}(p_s){\gamma }_{\mu }{P}_{L}u(p_b),$$

$${\mathcal{V}}_{sb}(q^2,M_D^2,m_I^2,m_R^2)$$

$$-2(C_{00}^{(IR)}{|}_{\mathrm{finite}}+C_{00}^{(RI)}{|}_{\mathrm{finite}}).$$

$$iA$$

Full text

1 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\to s\ell\ell$ transition has been drawing much interest for the last several years because of anomalies in $B\to K^{(*)}\mu\mu$ and $B_{s}\to\phi\mu\mu$ decays.

In particular SM predictions on the ratio of branching fractions

[TABLE]

are close to unity, signifying the lepton flavor universality (LFU) in the SM. However, the measurements at the LHCb for $K$ Aaij:2014ora and $K^{*}$ Aaij:2017vbb are lower than unity at $2.3-2.6\sigma$ level. Because the ratio (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\to K^{*}\mu^{+}\mu^{-}$ and branching fraction of $B_{s}\to\phi\mu^{+}\mu^{-}$ and $\Lambda_{b}\to\Lambda\mu^{+}\mu^{-}$, the deviations from the SM predictions increase as large as about 5$\sigma$ Altmannshofer:2017fio ; Capdevila:2017bsm ; Ciuchini:2017mik ; Alok:2017sui , which we will call $b\to s\mu\mu$ anomaly. At $m_{b}$ scale the $b\to s\ell\ell$ transition is described by the effective weak Hamiltonian

[TABLE]

where the relevant effective operators are

[TABLE]

In the SM the un-primed operators dominate the chirality-flipped primed ones. In the SM we obtain $C_{7\gamma}^{\rm SM}\simeq-0.294$, $C_{9}^{\rm SM}\simeq 4.20$, $C_{10}^{\rm SM}\simeq-4.01$ at $m_{b}$ scale Bobeth:1999mk ; Mahmoudi:2018qsk . The results from global fitting analyses Altmannshofer:2017fio ; Capdevila:2017bsm ; Ciuchini:2017mik ; Alok:2017sui show that sizable NP contributions to $C_{9}^{\mu}$ and/or $C_{10}^{\mu}$ can explain the $b\to s\mu\mu$ anomaly.

In this paper we consider a NP model with $C_{10}^{\mu,{\rm NP}}=0$, in which case the best fit value for $C_{9}^{\mu,{\rm NP}}$ is Capdevila:2017bsm

[TABLE]

with a SM pull of $5.8\sigma$. In addition to the SM gauge groups we introduce a new gauge symmetry $U(1)_{L_{\mu}-L_{\tau}}$ 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_{\mu}-L_{\tau}}$ gauge boson $Z^{\prime}$ couples to muon, it can make a contribution to the anomalous magnetic moment of muon $(g-2)_{\mu}$ Baek:2001kca ; Banerjee:2018eaf . But the $Z^{\prime}$ should be very light ($\lesssim 400$ MeV) to fully accommodate the discrepancy between the experiments and the SM predictions in the $(g-2)_{\mu}$ Altmannshofer:2014pba .

The model can also be extended to accommodate neutrino data Baek:2015mna ; Baek:2015fea ; Singirala:2018mio ; Asai:2018ocx . In Ref. Baek:2017sew we introduced a fermion dark matter (DM) model whose stability is originated from $U(1)_{L_{\mu}-L_{\tau}}$ symmetry Baek:2008nz . The model can also explain $b\to s\mu\mu$ 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 Crivellin:2015mga ; Belanger:2015nma ; Allanach:2015gkd ; Ko:2017quv ; Ko:2017yrd ; Ko:2017lzd ; Arnan:2016cpy ; Altmannshofer:2016jzy ; Kawamura:2017ecz ; Assad:2017iib ; Baek:2018aru ; Darme:2018hqg ; Barman:2018jhz ; Rocha-Moran:2018jzu ; Faisel:2018bvs ; Vicente:2018frk . In this paper we consider a “spin-flipped” version of the model in Ref. Baek:2017sew . We introduce two complex scalar fields $S$ and $X$: $S$ breaks the $U(1)_{L_{\mu}-L_{\tau}}$ symmetry spontaneously by developing vacuum expectation value (VEV) $\langle S\rangle$, while the lighter component $X_{I}$ is stable by the remnant discrete $Z_{2}$ symmetry after $U(1)_{L_{\mu}-L_{\tau}}$ symmetry is broken spontaneously and become a DM candidate. To explain the $b\to s\mu\mu$ anomaly as well in this model, we introduce a vector-like quark $Q$ which can mediate quark couplings to $Z^{\prime}$ boson. We will study the solution of the $b\to s\mu\mu$ 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\to s\mu\mu$, and consider low energy constraints including $C_{9}^{\ell,{\rm NP}}$, $\Delta m_{s}$ in $B_{s}-\bar{B}_{s}$ mixing, $B\to K^{(*)}\nu\bar{\nu}$, $b\to s\gamma$, the anomalous magnetic moment of muon $a_{\mu}$, and the loop-induced effective $Zb\bar{b}$ coupling. In Section 4 we consider dark matter phenomenology. Finally we conclude in Section 6. Loop functions are collected in Appendix A.

2 The model

We introduce a scalar dark matter candidate $X$ and a scalar boson $S$ which gives a mass to $U(1)_{L_{\mu}-L_{\tau}}$ gauge boson $Z^{\prime}$ after the symmetry is broken down spontaneously by the VEV of $S$. To couple the $Z^{\prime}$ gauge boson to the quarks we also introduce a vector-like $SU(2)_{L}$-doublet fermion $Q\equiv(U,D)^{T}$. Their charges under the $U(1)_{L_{\mu}-L_{\tau}}$ as well as those under the SM gauge groups are shown in Table 1.

The Lagrangian respecting the gauge symmetry and charge assignments in Table 1 is written as

[TABLE]

where $D_{\mu}$ is the covariant derivative, $i(=1,\cdots,3)$ is the quark-generation index, and $A_{\mu\nu}=\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}$ ($A=Z^{\prime},B$) is the field strength tensor. The scalar potential is in the form,

[TABLE]

The trilinear $\mu$ term allows a remnant discrete $Z_{2}$ symmetry after $S$ gets VEV and $U(1)_{L_{\mu}-L_{\tau}}$ is spontaneously broken. This local $Z_{2}$ symmetry Krauss:1988zc stabilises the DM candidate, which we assume the lighter component of $X$. The kinetic mixing angle $\chi$ is strongly constrained to a level of $O(10^{-3})$ by the DM direct search experiments Mambrini:2011dw . The non-vanishing $\chi$ does not help solving $b\to s\mu\mu$ anomaly because the SM gauge bosons allow only LFU couplings. In this paper we neglect this term for simplicity by setting $\chi\equiv 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_{\mu}-L_{\tau}}$ 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 $\mu$ term makes the complex scalar $X$ split into two real scalar fields $X_{R,I}$ defined by

[TABLE]

with masses-squared

[TABLE]

Assuming $\mu>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}$, $\lambda_{HX},\lambda_{SX}$ as free parameters, then we can write the parameters $m_{X}^{2}$ and $\mu$ in the Lagrangian as

[TABLE]

After $S$ gets VEV, the $U(1)_{L_{\mu}-L_{\tau}}$ gauge boson $Z^{\prime}$ obtains mass,

[TABLE]

where $g_{Z^{\prime}}$ is the gauge coupling constant of the $U(1)_{L_{\mu}-L_{\tau}}$ 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 Belanger:2015nma ; Altmannshofer:2014cfa .

In the unitary gauge we decompose the SM Higgs $H$ and the dark scalar $S$ as

[TABLE]

The stationary condition at the vacuum gives conditions

[TABLE]

Using the above conditions it is straightforward to obtain the scalar mass-squared matrix

[TABLE]

in the basis $(h,s)$. It is diagonalised by introducing mixing angle $\alpha_{H}$ to get the scalar mass eigenstates $(H_{1},H_{2})$

[TABLE]

where $H_{1}$ is identified with the SM-like Higgs boson with mass $m_{H_{1}}=125$ GeV. We will take $m_{H_{1}}$, $m_{H_{2}}$, and $\alpha_{H}$ as input parameters. Then the parameters $\lambda_{H,S}$ and $\lambda_{HS}$ are derived from them,

[TABLE]

where $c_{H}(s_{H})$ is an abbreviation of $\cos\alpha_{H}(\sin\alpha_{H})$. We require $\lambda_{HS}>-2\sqrt{\lambda_{H}\lambda_{S}}$ to stabilise the scalar potential at the electroweak (EW) scale.

In (5) 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 up-quark sector, i.e. $d_{iL}=d^{\prime}_{iL},u_{iL}=\sum_{j}V_{ji}^{*}u^{\prime}_{jL}$ where primes represent the mass eigenstates. Then the Yukawa interactions of the type $q-Q-X$ can be written in the form

[TABLE]

where $\lambda_{u_{i}}\equiv\sum_{j}V_{ij}\lambda_{j}$ and $\lambda_{d_{i}}\equiv\lambda_{i}$. We will also use the notation $\lambda_{d,s,b}$ for $\lambda_{d_{i}}$ and $\lambda_{u,c,t}$ for $\lambda_{u_{i}}(i=1,2,3)$. We will simply set $\lambda_{1}=0$ to remove the constraints related to the first generation quarks. Even in this case we see

[TABLE]

is induced. The induced $\lambda_{u}$ can generate NP contribution to $D^{0}-\bar{D}^{0}$ mixing. However, due to Cabibbo-suppressed contribution to $D^{0}-\bar{D}^{0}$ at least by ${\cal O}(\lambda^{2})$ where $\lambda(\approx 0.23)$ is the Cabibbo angle, the constraint from $D^{0}-\bar{D}^{0}$ can be always satisfied once the constraint from $B_{s}-\bar{B}_{s}$ is imposed Arnan:2016cpy . And we do not consider this constraint further.

The DM interacts with the SM fields through the Higgs-portal Lagrangian

[TABLE]

In this paper we will set $\alpha_{H}=0$ and $\lambda_{HX}=0$ to suppress the stringent constraint from the dark matter direct detection experiments via this Higgs portal interaction Baek:2011aa ; Baek:2012se .

3 $R_{K}^{(*)}$ and constraints

The $b\to s\mu^{+}\mu^{-}$ transition operator $O_{9}^{\mu}$ which can explain the $R_{K^{(*)}}$ anomaly is generated via the penguin diagrams shown in Fig. 1. The arrows represent color or lepton number flow. First, we calculate the one-loop effective $s(p_{s})-b(p_{b})-Z^{\prime}(q)$ ($q=p_{s}-p_{b}$) vertex. Assuming $m_{I},m_{R}$, and $M_{D}$ are at the EW scale ($\equiv M_{\rm EW}$), we can neglect terms proportional to external quark mass squareds, $m_{s(b)}^{2}/M_{\rm EW}^{2}(\ll 1)$. In this approximation, it is straightforward to get the effective vertex for diagrams (a) and (b) in Fig. 1:

[TABLE]

where we take the dimension of space-time integration $d$ to be $d\equiv 4-2\epsilon$ for positive infinitesimal $\epsilon$. The $C$’s are abbreviations for one-loop three-point functions defined in Hahn:1998yk ,

[TABLE]

where $k=0,12,00$ and $i,j=I,R(i\centernot=j)$. We will set $m_{s}=m_{b}\equiv 0$ in the calculation of the $C$-functions to be consistent with our approximation $m_{s(b)}^{2}/M_{\rm EW}^{2}\ll 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

[TABLE]

where $C_{00}|_{\rm finite}$ is the remaining finite part. Using the relation between the $U(1)_{L_{\mu}-L_{\tau}}$ charges, $q_{Q}=-q_{X}$, we can show that the sum of the two one-loop effective vertices is finite and given by

[TABLE]

where

[TABLE]

Now we can attach the external muon line in Fig. 1 to the $Z^{\prime}$ to get $C_{9}^{\mu}$. The full amplitude for $b\to s\mu\mu$ transition in Fig. 1 is given by

[TABLE]

where $p_{3}(p_{4})$ is outgoing (incoming) muon four-momentum. The term proportional to $q_{\mu}q_{\nu}$ vanishes because $\bar{u}(p_{3})\centernot qu(p_{4})=0$. Since $q\sim{\cal O}(m_{b})$ at most, we can set $q^{2}\equiv 0$ in the denominator of $Z^{\prime}$-propagator. In this case the effective vertex can be written in a simple analytic form:

[TABLE]

where $x_{I(R)}=m_{I(R)}^{2}/M_{D}^{2}$ and the loop function $k$ is defined in Appendix A. We note ${\cal V}_{sb}\to 0$, when $x_{R}\to x_{I}$. This can be understood as follows: in the limit $m_{R}\to m_{I}$, the two real scalars $X_{I}$ and $X_{R}$ merge into the original complex scalar $X$ as can be seen from (7). In this limit, the subset of the full Lagrangian which contributes the effective vertex $iV_{\rm eff}^{\mu}(q^{2})$,

[TABLE]

is invariant under local $U(1)_{L_{\mu}-L_{\tau}}$ symmetry. Note that the $Z^{\prime}$ mass term which breaks the $U(1)_{L_{\mu}-L_{\tau}}$ is not included in $\Delta{\cal L}$. Then the Ward-Takahashi identity dictates $V_{\rm eff}(q^{2}=0)=0$111This result holds also in case we keep the quark masses $m_{s}$ and $m_{b}$ because they respect the $U(1)_{L_{\mu}-L_{\tau}}$ symmetry. due to the absence of the tree-level $Z^{\prime}$-exchanging FCNC, leading to $V_{\rm eff}(q^{2})\propto q^{2}$ to all orders of perturbation theory  Peskin:1995ev . Since $V_{\rm eff}(q^{2}=0)=0$, we obtain ${\cal V}_{sb}=0$ in the limits $q^{2}\to 0$ and $m_{R}\to m_{I}$.

Now it is straightforward to get

[TABLE]

where the prime on the $k$ functions denotes a derivative with respect to the argument and we fixed the $U(1)_{L_{\mu}-L_{\tau}}$ charge of $\mu^{-}$ to be $+1$. A sizable mass splitting between $m_{R}$ and $m_{I}$ is favoured to generate $C_{9}^{\mu,{\rm NP}}$ which can explain the $b\to s\mu\mu$ anomaly. As a benchmark point in the parameter space, we choose $q_{X}=2$, $\alpha_{Z^{\prime}}=0.1$, $m_{Z^{\prime}}=700$ GeV, $\lambda_{s}\lambda_{b}^{*}=0.2$, $M_{D}=1$ TeV, $m_{I}=900$ GeV, and $m_{R}=3$ TeV, for which we get

[TABLE]

which is close to the best fit value in (4) to solve the $b\to s\mu\mu$ anomaly. Fig. 2(a) shows $C_{9}^{\mu,{\rm NP}}$ as a function of $m_{I}$ for three different values of $m_{R}=1,2,3$ TeV (from above) with $m_{Z^{\prime}}=700$ GeV, $M_{D}=1$ TeV, $q_{X}=2$, $\alpha_{Z^{\prime}}=0.1$, and $\lambda_{s}\lambda_{b}^{*}=2$. As can be seen from (35), $|C_{9}^{\mu,{\rm NP}}|$ becomes larger as the mass splitting $m_{R}-m_{I}$ increases. Fig. 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 Fig. 2(a). We can see that the effect of the vector-like $D$-quark decouples as $M_{D}$ increases. In both cases, smaller $m_{I}$ is favored to obtain larger $|C_{9}^{\mu,{\rm NP}}|$.

Photon- or $Z$-penguin diagrams similar to $Z^{\prime}$-penguin diagrams in Fig. 1 but with $Z^{\prime}$ replaced by photon or $Z$-boson can contribute to $C_{9}^{\mu,{\rm NP}}$. Their couplings to leptons are flavour-universal and they also contribute to $C_{9}^{e,{\rm NP}}$ and $C_{9}^{\tau,{\rm NP}}$ 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_{\rm eff}(q^{2})\propto q^{2}$ above. Here the conserved $U(1)$ symmetries are the $U(1)$-electromagnetism, $U(1)_{\rm 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_{\mu}-L_{\tau}}$ is conserved or not, the argument applies even when $m_{I}\centernot=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 $q^{2}/M_{Z}^{2}$ with $q^{2}\sim m_{b}^{2}$. We obtain the photon penguin contribution to be

[TABLE]

where $x_{I(R)}=m_{I(R)}^{2}/M_{D}^{2}$ and the loop function $Q_{1}(x)$ is listed in (75). For the benchmark point $\lambda_{s}\lambda_{b}^{*}=0.2$, $M_{D}=1$ TeV, $m_{I}=900$ GeV, and $m_{R}=3$ TeV, we get

[TABLE]

which is about three orders of magnitude smaller than the $Z^{\prime}$ contribution to $C_{9}^{\mu,{\rm NP}}$ in (36). 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 $|\lambda_{s}\lambda_{b}^{*}|$ is the most strongly constrained by the measurements of the mass difference $\Delta m_{s}$ for $B_{s}-\bar{B}_{s}$ mixing. Fig. 3 shows one-loop box diagrams for $B_{s}-\bar{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,

[TABLE]

The Wilson coefficient $C_{1}$ can be decomposed into the SM and the NP contributions

[TABLE]

The SM contribution at the electroweak scale is obtained by box diagrams with $W$-boson and $t$-quark running inside the loop:

[TABLE]

where $x_{t}=m_{t}^{2}/m_{W}^{2}\approx 4.64$ and the loop function $S_{0}(x_{t})$ can be found in Buchalla:1995vs . The NP diagrams shown in Fig. 3 give

[TABLE]

where $x_{i}=m_{i}^{2}/M^{2}_{D},\,(i=R,I)$. We note that this result is non-vanishing, different from a single real DM contribution which vanishes Arnan:2016cpy . 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}-\bar{B}_{s}$ system gives a constraint on the value of $C_{1}^{\rm NP}$:

[TABLE]

at 2$\sigma$ confidence level Arnan:2016cpy . For the benchmark point $m_{I}=900$ GeV, $M_{D}=1$ TeV and $m_{R}=3$ TeV, we get

[TABLE]

which is about an order of magnitude smaller than the SM prediction. This point satisfies the constraint (43).

Fig. 4 shows a contour plot of $C_{1}^{\rm NP}$ with contour lines $0.03,0.1,0.3$, and $0.6$ in the unit of $10^{-11}\,{\rm GeV^{-2}}$ in the $(m_{I},M_{D})$-plane for $\lambda_{s}\lambda_{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. The green and yellow region is excluded by (43). The plot shows that $C_{1}^{\rm NP}$ is always positive in our model and the constraint (44) is easily satisfied when new particles are at TeV scale.

Another possible constraint on the model parameters comes from the experimental measurements of the inclusive branching fraction of radiative $B$-decay, $\bar{B}\to X_{s}\gamma$ Amhis:2016xyh ,

[TABLE]

For this process the SM prediction has been calculated up to NNLO QCD corrections Misiak:2015xwa , which predict,

[TABLE]

The NP contribution to ${\cal B}(\bar{B}\to X_{s}\gamma)$ can be obtained by calculating the Wilson coefficients $C_{7\gamma,8g}$ from the diagrams in Fig. 5:

[TABLE]

where $e_{D}=-1/3$ is the electric charge of the vector-like down-type quark $D$ and $x_{i}=m_{i}^{2}/M_{D}^{2}\,(i=I,R)$. The loop-fuction $J_{1}(x)$ is listed in the Appendix A. From the prediction including NP contribution to $C_{7\gamma(8g)}$ Misiak:2015xwa , (45) and (46), we obtain the constraint

[TABLE]

at 2$\sigma$ level. For the benchmark point $m_{I}=900$ GeV, $M_{D}=1$ TeV, and $m_{R}=3$ TeV, we obtain

[TABLE]

which is about two orders of magnitude less than the current bound (48). Fig. 6 shows a contour plot for the combination $C_{7\gamma}+0.24\,C_{8g}$ with contour lines $-10^{-5},-5\times 10^{-5}$, and $-10^{-4}$ in the $(m_{I},M_{D})$-plane for $\lambda_{s}\lambda_{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\gamma}+0.24\,C_{8g}$ is less sensitive to $m_{R}$ than $C_{1}^{\rm NP}$ of $B_{s}-\bar{B}_{s}$ is. The entire region considered is allowed by ${\cal B}^{\rm exp}(\bar{B}\to X_{s}\gamma)$.

The NP diagrams for semi-leptonic decay $B\to K^{(*)}\nu\bar{\nu}$ is obtained when the external muon lines are replaced with neutrino lines in Fig. 1. The effective Hamiltonian for the decay is

[TABLE]

where

[TABLE]

We obtain the non-vanishing coefficients

[TABLE]

We note that the diagram with $Z^{\prime}$ replaced by $Z$ vanishes in the $q^{2}\to 0$ limit, showing that the $Z^{\prime}$ contribution is dominant. The current experimental bounds on the ratios of branching fractions

[TABLE]

are

[TABLE]

In our model these ratios are predicted to be

[TABLE]

where we used $C_{L}^{33,{\rm NP}}=-C_{L}^{22,{\rm NP}}\centernot=0$ while all the other components are vanishing. We can see the interference terms cancel each other out. Considering $C_{L}^{22,{\rm NP}}=C_{9}^{\mu,{\rm NP}}/2\approx-0.6$ to explain the $b\to s\mu\mu$ anomaly and $C_{L}^{\rm SM}\approx-6.35$, we predict

[TABLE]

showing the deviation from the SM is very small partly due to the cancellation of the interference terms.

The gauged $U(1)_{L_{\mu}-L_{\tau}}$ model is well known to generate the sizable anomalous magnetic moment of muon, $a_{\mu}=(g-2)_{\mu}/2$ via the $Z^{\prime}$-exchanging one-loop diagram Baek:2001kca . The $Z^{\prime}$ contribution can explain the long-standing discrepancy between the experimental measurements Bennett:2006fi and the SM predictions Kurz:2014wya :

[TABLE]

The effective Hamiltonian for $a_{\mu}$ is

[TABLE]

The NP contribution at one-loop level is calculated to be

[TABLE]

which in the limit, $m_{\mu}^{2}\ll m_{Z^{\prime}}^{2}$, approximates

[TABLE]

For the benchmark point $m_{Z^{\prime}}=700$ GeV, $\alpha_{Z^{\prime}}=0.1$, we get

[TABLE]

which is consistent with (57) within 3$\sigma$. In the minimal $U(1)_{L_{\mu}-L_{\tau}}$ model the region in the $(m_{Z^{\prime}},\alpha_{Z^{\prime}})$ plane which can explain $\Delta a_{\mu}$ at 2$\sigma$ level is excluded by the bound from the measurement of neutrino trident production, $\nu_{\mu}N\to\nu_{\mu}N\mu^{+}\mu^{-}$, when $m_{Z^{\prime}}\gtrsim 400$ MeV Altmannshofer:2014pba . Fig. 7 shows the constraint from neutrino trident production and muon $g-2$ in $(m_{Z^{\prime}},\alpha_{Z^{\prime}})$ plane. The grey region is disfavoured by the neutrino trident production experiments at 2$\sigma$. The region between the two green lines is favoured by the current discrepancy $\Delta a_{\mu}$ at 2$\sigma$, but it is excluded by the neutrino trident production experiments. The red “X” mark represents the benchmark point $m_{Z^{\prime}}=700$ GeV and $\alpha_{Z^{\prime}}=0.1$.

The new particles in the model also generates one-loop effective $Zf\bar{f}$-vertex $(f=s,b)$. Since $Zb\bar{b}$ vertex has been more precisely determined by the LEP experiment, we consider the constraint only from $Zb\overline{b}$. The $Zb\bar{b}$ vertex is written in the form,

[TABLE]

where tree-leve values for the couplings are $g_{L}^{b}({\rm tree})=-1/2-e_{d}{s_{w}}^{2}$, $g_{R}^{b}({\rm tree})=-e_{d}{s_{w}}^{2}$, $e_{d}=-1/3$, and ${s_{w}}^{2}\approx 0.23$. The deviation of $g_{R(L)}^{b}$ from the SM prediction obtained from a global fit is Ciuchini:2013pca 222We choose more conservative result in Ciuchini:2013pca .

[TABLE]

with a correlation coefficient of $+0.9$. In our model the NP contributions to $g_{L}^{b}$ is obtained to be333The new particles couple to only $b_{L}$ and do not generate $g_{R}^{b}$.,

[TABLE]

where $x_{I(R)}=m_{I(R)}^{2}/M_{D}^{2}$ and we can take $q^{2}=m_{Z}^{2}$. The loop function $Q_{1}(x)$ is listed in (75). We notice that the loop function is the same with the one for the photon penguin diagram of $b\to s\ell\ell$ in (37). 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 (37). So they should be proportional to each other. The above result can be compared with (63). In Fig. 8 we show error ellipses at 1, 2, and 3$\sigma$ confidence level in $(\delta g_{R}^{b},\delta g_{L}^{b})$ plane. The vertical red line segment is obtained by randomly scanning $g_{L}^{b,{\rm NP}}$ in the ranges

[TABLE]

The model predicts $\delta g_{L}^{b,{\rm NP}}$ in the range $(-4.9\times 10^{-4},0)$ and satisfies (63) at 3$\sigma$ 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\to D\bar{D}(U\bar{U})$ at the LHC if their masses are within the LHC reach. These production processes are similar to those considered in Baek:2016lnv ; Baek:2017ykw where they were analysed in detail. Roughly $M_{D(U)}\lesssim 1$ TeV are excluded. And we impose $M_{D}\geq 1$ TeV.

4 The dark matter

In this section we identify the main channel and the favoured parameter region to give the observed DM relic density, $\Omega_{\rm DM}h^{2}=0.1199\pm 0.0022$ Ade:2015xua . 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\sim 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\sim m_{I}/25$ Kolb:1990vq . Then the relic density is roughly related with the annihilation cross section at freeze-out temperature as

[TABLE]

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 Cline:2013gha . The scalar potential in the minimal model has terms

[TABLE]

The DM mass $m_{D}$ is obtained by $m_{D}^{2}=\mu_{D}^{2}+\lambda_{hD}v_{H}^{2}/2$. The DM annihilation occurs through $DD\to h\to\mathrm{SM}\,\mathrm{SM}$ or $DD\to hh$. Both processes are controlled by the Higgs portal coupling $\lambda_{hD}$, which is strongly restricted by the direct detection experiments  Akerib:2016vxi ; Cui:2017nnn ; Aprile:2018dbl . As a consequence the model is strongly constrained, ruling out $m_{D}\lesssim 1$ TeV region as a single-component DM Athron:2017kgt .

In our model, however, there are many model parameters involved in the DM-Higgs couplings as can be seen in (26), which allows the direct detection constraint onthe Higgs portal interaction to be lifted by setting $\alpha_{H}=\lambda_{HX}=0$ to remove $H_{1}X_{I}^{2}$ 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}\to Z^{\prime}Z^{\prime}$ and $X_{I}X_{I}\to q\bar{q}$. Barring the Higgs portal $X_{I}$ interaction with $H_{1}$, they are dominant processes. In Fig. 9 we show representative diagrams for the two annihilation channels. We implemented our model to the micrOMEGAs package Belanger:2013oya to evaluate the DM relic density and direct detection cross section.

We first consider the scenario in which the diagrams of type (a) in Fig. 9 play a major role. The process $X_{I}X_{I}\to Z^{\prime}Z^{\prime}$ dominates the DM annihilations as long as it is kinematically open, $\alpha_{Z^{\prime}}$ is not too small ($\alpha_{Z^{\prime}}\gtrsim 10^{-6}$), and $m_{H_{2}}$ is not much larger than TeV scale ($m_{H_{2}}\lesssim 5$ TeV). In this case we obtain typically $\langle\sigma v(X_{I}X_{I}\to Z^{\prime}Z^{\prime})\rangle\gg\langle\sigma v(X_{I}X_{I}\to q\bar{q}^{\prime})\rangle$ for the thermal-averaged annihilation cross sections. Given that we set $\alpha_{H}=\lambda_{HX}=0$, the process $X_{I}X_{I}\to Z^{\prime}Z^{\prime}$ is controlled by the dark Higgs interaction and the dark gauge interaction.

The former interaction is given by $\lambda_{S}v_{S}-\sqrt{2}\mu$, and the latter by $g_{Z^{\prime}}v_{S}$. Both are sensitive to $C_{9}^{\mu,{\rm NP}}$ in (35). Fig. 10 shows contour lines of $\Omega_{\rm DM}h^{2}=0.1199$ in $(m_{I},\alpha_{Z^{\prime}})$ (left panel) and $(m_{I},m_{R})$ (right panel) plane for $m_{Z^{\prime}}=10,100,700$ GeV. For the other parameters we take the benchmark values: $q_{X}=2$, $\lambda_{s}\lambda_{b}^{*}=0.2$, $M_{D}=1$ TeV, and $m_{H_{2}}=2$ TeV. We set $m_{R}=3$ TeV for the left panel and $\alpha_{Z^{\prime}}=0.1$ for the right panel. We also fixed $\alpha_{H}=0$, $\lambda_{X}=1$, $\lambda_{HX}=0$, and $\lambda_{SX}=0.2$. At this stage we do not impose constraints other than $m_{R}>m_{I}$. Numerically we have checked that $\langle\sigma v(X_{I}X_{I}\to Z^{\prime}Z^{\prime})\rangle$ is much larger than $\langle\sigma v(X_{I}X_{I}\to q\bar{q}^{\prime})\rangle$ for the points in Fig. 10.

The process $X_{I}X_{I}\to Z^{\prime}Z^{\prime}$ can occur in the early universe even when $m_{I}<m_{Z^{\prime}}$ if the mass difference is not too large. This can be seen in the steep lines corresponding to $m_{I}<m_{Z^{\prime}}$ in the left panel of Fig. 10. This is possible when the DMs move fast and their center of mass energy exceeds twice the $m_{Z^{\prime}}$: $\sqrt{s}>2m_{Z^{\prime}}$. To produce on-shell $Z^{\prime}$-pair, the relative velocity of DM pair in the CM-frame should satisfy

[TABLE]

For example, for $m_{I}=90$ GeV and $m_{Z^{\prime}}=100$ GeV, we obtain $v\geq 0.87$. The DM should be quite relativistic and the thermally averaged annhilation cross section is Boltzmann-suppressed. When $X_{I}X_{I}\to Z^{\prime}Z^{\prime}$ is kinematically open for non-relativistic $X_{I}$, the process is sensitive to the dark gauge coupling $\alpha_{Z^{\prime}}$. 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 $\mu-$term in (26). This shows that the process is also sensitive to the mass-squared difference, $m_{R}^{2}-m_{I}^{2}$, by (9). When $m_{I}$ is not close to the resonance region, the $s$-wave annihilation cross section for the $X_{I}X_{I}\to H_{2}\to Z^{\prime}Z^{\prime}$ channel is in the form

[TABLE]

When $\lambda_{SX}$ is not too large, the larger mass squared difference $m_{R}^{2}-m_{I}^{2}$ and the smaller $m_{Z^{\prime}}^{2}$, i.e. the larger $\mu$, the larger $\sigma v$ is obtained.

As $m_{I}$ approaches 1 TeV, it is close to the resonance region $m_{H_{2}}\approx 2m_{I}$ and the cross section increases rapidly, virtually independent of $\alpha_{Z^{\prime}}$. This explains almost vertical parts of the curves near $m_{I}=1$ TeV.

In the right panel of Fig. 10, we can see that the regions $m_{R}\approx m_{I}$ also give the correct relic density. This occurs due to the coannihilation processes $X_{R}X_{I}\to Z^{\prime}H_{i}(q\bar{q},\ell\bar{\ell})$ Baek:2017sew . As we saw in (35), the NP contribution to $C_{9}^{\mu,{\rm NP}}$ is suppressed. And the coannihilation mechanism for the DM relic density is not favoured as a solution to $b\to s\mu\mu$ anomaly. This shows a strong interplay between the flavour physics and DM phenomena Baek:2001nz ; Baek:2004et ; Baek:2005wi ; Baek:2005di ; Baek:2015fma ; Baek:2016kud ; Baek:2017sew ; Baek:2017ykw , which will be discussed in the next Section in more detail. Once $X_{I}X_{I}\to Z^{\prime}Z^{\prime}$ is kinematically open near $m_{I}=m_{Z^{\prime}}$, it dominates the annihilation processes, which is not so sensitive to $m_{R}$.

Now we consider the parameter space where type (b) diagrams in Fig. 9 dominate the annihilation cross section. Fig. 11 shows the dependence of the dark Yukawa coupling $\lambda_{b}$ as a function of $m_{I}$ to give constant $\Omega_{\rm 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^{\prime}}=10$ TeV and $m_{H_{2}}=10$ TeV to suppress the $X_{I}X_{I}\to Z^{\prime}Z^{\prime}$ channel. We set $\lambda_{s}=0.4$. For the other fixed parameters, we take the same values with those for Fig. 10. The $\lambda_{b}$ required to give $\Omega_{\rm DM}h^{2}$ changes sharply near $m_{t}/2$ and $m_{t}$. This occurs due to the processes $X_{I}X_{I}\to c\bar{t}(\bar{c}t)$ and $X_{I}X_{I}\to t\bar{t}$, respectively. Their $s$-wave annihilation cross sections are given by

[TABLE]

where we have neglected the mass of charm quark. Note that both are proportional to $m_{t}^{2}$. The contribution from $X_{I}X_{I}\to b\bar{b}$, being proportional to $m_{b}^{2}$, is negligible compared to the above two processes444In the limit $m_{b}\to 0$, the leading term in $v$ is proportional to $v^{4}$, i.e. $d$-wave.. When $m_{I}\approx M_{D}$, the near vertical lines are due to the coannihilation processes such as $D\bar{D}(U\bar{U})\to gg,q\bar{q},Zg$ $(q=u,d,s,c,b,t)$, $U\bar{D}\to W^{+}g$, and $X_{I}D\to sg$. They are sensitive to the SM $SU(3)_{C}$ gauge coupling and independent of $\lambda_{b}$. If we require $\lambda_{b}$ to be of order 1, to give the correct relic density $M_{D}$ should not be much heavier than ${\cal O}(1)$ TeV, increasing the prospect of producing $D$ or $U$ at the LHC.

5 Interplay between the $b\to s\mu\mu$ anomaly and the dark matter

Now we investigate whether the parameter space which solves the $b\to s\mu\mu$ 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 Akerib:2016vxi , PANDA Cui:2017nnn , and XENON1T Aprile:2018dbl .

Fig. 12 shows plots for $C_{9}^{\mu,{\rm NP}}$ which solves the $b\to s\mu\mu$ anomaly at 1$\sigma$ (dark blue) and 2$\sigma$ (light blue) in $(m_{I},\alpha_{Z^{\prime}})$ plane. We take $m_{Z^{\prime}}=700,100,10$ GeV (from the left panel). We fixed $q_{X}=2$, $\lambda_{s}\lambda_{b}^{*}=0.2$, $M_{D}=1$ TeV, and $m_{R}=3$ TeV. They are superimposed with the constant lines for $\Omega_{\rm DM}h^{2}=0.1199$ shown in Fig. 10. The grey regions are excluded by the neutrino trident production experiments at 2$\sigma$ level. We checked that the other low energy experiments do not further constrain the allows regions for the $C_{9}^{\mu,{\rm NP}}$ and the relic density. Neither does the direct detection experiments affect the plots in Fig. 12 because i) the $Z^{\prime}$ couples to the quarks at one-loop level and ii) more importantly only inelastic upward scattering $X_{I}q\to X_{R}q$ can occur for $Z^{\prime}$ interaction, which is forbidden kinematically. We can see that the $b\to s\mu\mu$ anomaly can be resolved at 1$\sigma$ for $m_{Z^{\prime}}=700$ GeV and the current dark matter can be accommodated at the same time. For smaller $Z^{\prime}$ masses, $m_{Z^{\prime}}=100,10$ GeV, $C_{9}^{\mu,{\rm NP}}$ becomes too large and $b\to s\mu\mu$ anomaly can be explained only at 2$\sigma$ 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 Fig. 13 we show the $C_{9}^{\mu,{\rm NP}}$ for $b\to s\mu\mu$ 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^{\prime}}=700,100,10$ GeV. As in the previous case we fixed $q_{X}=2$, $\alpha_{Z^{\prime}}=0.1$, $\lambda_{s}\lambda_{b}^{*}=0.2$, and $M_{D}=1$ TeV. The other fixed parameters are the same with those for Fig. 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}-\bar{B}_{s}$ mixing, neutrino trident production, or direct detection of dark matter. For each case there is intersection region of the required $C_{9}^{\mu,{\rm NP}}$ and the correct relic density free from other experimental constraints. The region occurs near the kinematic threshold of $X_{I}X_{I}\to Z^{\prime}Z^{\prime}$. Relatively large $m_{R}$ compared to $m_{I}$ is also required to get sizable $C_{9}^{\mu,{\rm NP}}$.

Now we consider the impact of $\lambda$ couplings on the dark matter and the $C_{9}^{\mu,{\rm NP}}$. For this purpose we suppress the dark gauge contributions to them by decoupling $Z^{\prime}$ as in Fig. 11. To decouple we assume $Z^{\prime}$ is heavy: for our purpose it is enough to set $m_{Z^{\prime}}=2$ TeV as in Fig. 11. Fig. 14 shows the results of this setting in $(m_{I},\lambda_{b})$ plane. We fixed $q_{X}=2$, $\alpha_{Z^{\prime}}=0.1$, $\lambda_{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 Fig. 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 $\Delta m_{s}$ of $B_{s}-\bar{B}_{s}$ system. We can see that the direct detection experiments and the $B_{s}-\bar{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\to s\mu\mu$ anomaly can be explained simultaneously for the electroweak scale DM. We notice that the contribution to $C_{9}^{\mu,{\rm NP}}$ is not easily decoupled for very heavy $m_{Z^{\prime}}$ and $m_{R}$ due to large mass splitting $m_{R}-m_{I}$. Eventually as $m_{Z^{\prime}}$ and/or $m_{R}$ becomes even heavier, their impact on $C_{9}^{\mu,{\rm NP}}$ 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.

6 Conclusions

We considered a new physics model with the $U(1)_{L_{\mu}-L_{\tau}}$ 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\to s\mu\mu$ 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$ whose VEV breaks the dark $U(1)_{L_{\mu}-L_{\tau}}$ gauge symmetry spontaneously. Since the field $Q$ is vector-like under the gauge group, the model is free of the gauge anomaly. Their charges are assigned in such a way that after the dark Higgs $S$ gets a VEV, there still remains a remnant discrete $Z_{2}$ symmetry of $U(1)_{L_{\mu}-L_{\tau}}$. After the $U(1)_{L_{\mu}-L_{\tau}}$ is broken down to the $Z_{2}$, the complex field $X$ is split into the two real scalar fields $X_{I}$ and $X_{R}$. The $X_{I}$, being the lightest $Z_{2}$-odd particle, is the dark matter candidate in the model. The other $Z_{2}$-odd particles are the $X_{R}$, and the $Q$.

We identified some benchmark points which can explain both the $b\to s\mu\mu$ anomaly and the relic abundance of dark matter. We checked that they avoid the known constraints from the $B$-meson decays, the $K$-meson decays, the LEP experiments, the measurements of the anomalous magnetic moment of muon, and also the experiments of the direct detection of dark matter. When the $X_{I}X_{I}\to Z^{\prime}Z^{\prime}$ annihilation diagrams dominate, the $Z^{\prime}$, $X_{R}$, $X_{I}$, $Q$ with the electroweak scale masses can explain both the DM relic density and the $b\to s\mu\mu$ anomaly while satisfying the current constraints. Since the $C_{9}^{\mu,{\rm NP}}$ for $b\to s\mu\mu$ anomaly requires a sizable mass difference $m_{R}-m_{I}$, the coannihilation region $m_{R}\approx m_{I}$ is ruled out in this case. When the $X_{I}X_{I}\to q\bar{q}$ mode dominates the relic density calculation, the $Z^{\prime}$, $X_{R}$, and the dark Higgs can be much heavier than the electroweak scale ($\gtrsim 10$ TeV) while the DM, $Q$ are still at the TeV scale. When $\lambda_{b}\sim{\cal O}(1)$ the $C_{9}^{\mu,{\rm NP}}$ is still sizable because it is not easily decoupled with large $m_{R}-m_{I}$. The DM relic density can be explained by either the $X_{I}X_{I}\to t\bar{t}$ channel or the $Q\bar{Q}\to{\rm SM}\;{\rm SM}$ channel.

The model is a spin-flipped version of the model considered in Baek:2017sew and shares some results in common. In both models the $Z^{\prime}$-penguin diagram can accommodate the required $C_{9}^{\mu,{\rm NP}}$ to explain the $b\to s\mu\mu$ anomaly and the dark matter candidate can explain the current relic density of the universe. The strongest flavour constraint comes from the mass difference in the $B_{s}-\bar{B}_{s}$ system in both models. Compared with Baek:2017sew , we included some new constraints such as the photon ($Z$-boson) penguin diagrams, the new physics contributions to the $Zb\bar{b}$ vertex, and extended the discussion on the dark matter.

Appendix A Loop functions

The loop function $k(x)$ is defined as

[TABLE]

We obtain $k(1)=1$. When the number of arguments is greater than one, the loop function is defined recursively by

[TABLE]

From this definition, for example, we get

[TABLE]

with $k(1,1)=3/2$. The loop function $J_{1}(x)$ for $b\to s\gamma$ is

[TABLE]

We have $J_{1}(1)=1/24$. The loop function for the photon-penguin of $b\to s\ell\bar{\ell}$ the effective $Zb\bar{b}$ vertex is

[TABLE]

We get $Q_{1}(1)=1/8$.

Note added: After finalizing the manuscript we received a paper considering $b\to s\mu\mu$ anomaly in a similar but a different setting Hutauruk:2019crc . Their model does not have $\mu$ 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 $\mu$ term.

Acknowledgements.

This work was supported by the National Research Foundation of Korea(NRF) grant funded by the Korea government(MSIT) (Grant No. NRF-2018R1A2A3075605).