Gauged $U(1)_{L_\mu -L_\tau }$ scotogenic model in light of $R_{K^{(*)}}$ anomaly and AMS-02 positron excess

Abstract

We study the gauged $U(1)_{L_\mu -L_\tau }$ scotogenic model with emphasis on latest measurement of LHCb $R_{K^{(*)}}$ anomaly and AMS-02 positron excess. In this model, neutrino masses are induced at one-loop level with $Z_2$-odd particles, i.e., right-handed neutrinos $N_\ell (\ell =e,\mu ,\tau )$ and inert scalar doublet $\eta$ inside the loop. Meanwhile, the gauged $U(1)_{L_\mu -L_\tau }$ symmetry is broken spontaneously by the scalar singlet S, resulting to the massive gauge boson $Z'$. Provided certain couplings to quarks induced by heavy vector-like quarks, the gauge boson $Z'$ would contribute to the transition $b\rightarrow s \mu ^+\mu ^-$, hence explain the $R_{K^{(*)}}$ anomaly. As for the Majorana fermion DM N, the gauge boson $Z'$ and the singlet Higgs $H_0$ will generate various annihilation channels, among which the $NN\rightarrow Z'Z'$ and $NN\rightarrow Z'H_0(\rightarrow Z'Z')$ channel could be used to interpret the AMS-02 positron excess. We give a comprehensive analysis on model parameter space with consider various current constraints. The combined analysis shows that the $R_{K^{(*)}}$ anomaly and AMS-02 positron excess can be explained simultaneously.

1 Introduction

Tiny neutrino mass and non-baryonic dark matter (DM) are the two missing pieces of standard model (SM). An appealing pathway to link them together is the scotogenic model [1, 2], which realizes tiny neutrino mass via radiative process [3,4,5] with DM running in the loop. Along this idea, various possibilities [6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38] have been proposed.¹ However, in the minimal version of scotogenic model, the parameter space of fermionic DM required by relic abundance is tightly constrained by lepton flavor violation (LFV) processes [40,41,42]. Such tension motivates the suggestions that extend the original model with a new U(1) gauge group. The contradiction is then relaxed due to new available annihilation channels via exchanging of gauge or Higgs boson associated with the U(1) group. Particularly, comparing with gauged $U(1)_{B-L}$ model [43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59], gauged $U(1)_{L_\mu -L_\tau }$ model [60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77]) has less stringent constraints due to the fact that corresponding gauge boson $Z^\prime$ does not couple to SM quarks and electron directly. It is worthy to note that a light $Z^\prime \sim {\mathcal {O}}(100)~\mathrm{MeV}$ with gauge coupling $g^\prime \sim 10^{-3}$ is suitable to interpret the muon $g-2$ anomaly [78, 79].

Table 1 The particle content and the charge assignment under $SU(2)_L \times U(1)_Y \times U(1)_{L_\mu -L_\tau }\! \times Z_2$

Except the evidences from neutrino mass and DM, the hint from flavor physics may also call for the physics beyond standard model (BSM). Recently, a tentative evidence indicates lepton flavor universality (LFU) violation has been reported by LHCb Collaboration in the semi-leptonic decays of the B meson. The latest result gives [80, 81]

$$\begin{aligned} R_K= & {} \frac{\text {Br}(B\rightarrow K \mu ^+\mu ^-)}{\text {Br}(B\rightarrow K e^+e^-)}\nonumber \\= & {} 0.846^{+0.060+0.016}_{-0.054-0.014}\quad \text {for}\quad 1.1<q^2<6~\mathrm{GeV}^2, \end{aligned}$$

(1)

where $q^2$ is squared momentum of the leptonic system. This result presents $2.5\sigma$ deviation with respect to the SM prediction [82]. In addition, observations of a tension in angular observable, such as $P_5^\prime$ in the decay $B\rightarrow K^*\mu ^+\mu ^-$ and angular distribution in the decay $B_s^0\rightarrow \phi \mu ^+\mu ^-$, have been announced by LHCb [83] and Belle [84] as well. Since the decay process $b\rightarrow s\ell \ell (\ell =e,\mu )$ is involved in the aforementioned B anomalies, a new physics contribution to the corresponding Wilson coefficient is able to explain such anomaly [85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103,104,105,106,107,108,109,110,111,112,113,114,115]. Especially, the LFU violation and angular anomaly in $b\rightarrow s\mu ^+\mu ^-$ indicate that the new physics particles may prefer couple to muon rather than to electron, which is actually an intrinsic feature of $U(1)_{L_\mu -L_\tau }$ gauge boson $Z^\prime$ [116,117,118,119,120,121,122,123,124,125,126,127].

In addition, the latest direct detection experiments, such as LUX [128], XENON1T [129, 130] and PandaX-II [131] remains for DM signal. Since these experiments based on DM-hadron interaction, the null results of DM direct detection signal may suggest that the DM-hadron interaction is at least suppressed or even better vanishing. On the other hand, the indirect detection experiments, such as PAMELA [132], Fermi-LAT [133], and AMS-02 [134,135,136], have reported a significance positron fraction excess in the cosmic-ray, while no obvious antiproton excess is observed [137]. This can be interpreted by DM dominantly annihilating into leptonic final states [138,139,140,141,142].² Therefore, the current direct and indirect detection experiments seem to advocate ‘leptophilic DM’ [151]. And $U(1)_{L_\mu -L_\tau }$ DM is clearly a good choice [62].

With holding such benefits, we thus consider the phenomenology of the gauged $U(1)_{L_\mu -L_\tau }$ scotogenic model [65] in this paper, which has been studied in Refs. [65, 66] with emphasizing neutrino mass and dark matter properties for light $Z^\prime$. According to Ref. [66], $50~\mathrm{MeV}\lesssim M_{Z'}\lesssim 400~\mathrm{MeV}$ with $3\times 10^{-4}\lesssim g'\lesssim 10^{-3}$ is required to explain $(g-2)_\mu$ and DM relic density. Lately, searches for $Z'$ in the $e^+e^-\rightarrow \mu ^+\mu ^- Z', Z'\rightarrow \mu ^+\mu ^-$ channel by BABAR has rejected $M_{Z'}>212~\mathrm{MeV}$ with $g'\gtrsim 7\times 10^{-4}$ [152]. Although not fully excluded at current, the future experiments, such as Belle-II and SPS, are able to exclude the whole low $Z'$ mass region favored by $(g-2)_\mu$ and DM [153,154,155,156]. Therefore, we move on to the high mass region $M_{Z'}\gtrsim 10~\mathrm{GeV}$ and specially focus on $R_{K^{(*)}}$ anomaly and AMS-02 positron excess, which has not been discussed in previous studies [65, 66].

This paper is organized as following. In Sect. 2, we briefly review the gauged $U(1)_{L_\mu -L_\tau }$ scotogenic model. Interpretation of the $R_{K^{(*)}}$ anomaly is discussed in Sect. 3, and relative constraints on $Z'$ in the high mass region are summarised in Sect. 3.1. A detail scanning on the DM parameter space under constraints from relic density and direct detection are performed in Sect. 4. In Sect. 5, we provide benchmarks to fit positron excess by using latest AMS-02 measurement and discuss relative constraints imposed by indirect detections. Finally, Sect. 6 is devoted to conclusion.

2 The model

2.1 Model setup

The gauged $U(1)_{L\mu -L_\tau }$ scotogenic model was proposed in Ref. [65]. Comparing with the original scotogenic model [2], this model further introduces one additional singlet scalar S charged $+1$ under $U(1)_{L_\mu -L_\tau }$ to break this symmetry spontaneously. The particle content and the charge assignment under $SU(2)_L \times U(1)_Y \times U(1)_{L_\mu -L_\tau } \times Z_2$ are shown in Table 1. Here, the right-handed neutrinos $N_{\ell }(\ell =e,\mu ,\tau )$ and inert scalar doublet $\eta$ have $Z_2$ odd charge, thus the lightest of them can be a dark matter candidate. In this paper, we consider fermion DM. Neutrino masses are generated at one-loop level via $N_\ell$ and $\eta$ propagating in the loop as the original scotogenic model.

The scalar potential involving scalar doublet $\Phi$, scalar singlet S and inert scalar doublet $\eta$ is [65]

$$\begin{aligned} V= & {} + \mu _\Phi ^2 |\Phi |^2 + \mu _\eta ^2 |\eta |^2 + \mu _S^2 |S|^2 + \frac{1}{2} \lambda _1 |\Phi |^4 + \frac{1}{2} \lambda _2 |\eta |^4 \nonumber \\&+ \lambda _3 |\Phi |^2 |\eta |^2 + \lambda _4 |\Phi ^\dagger \eta |^2 +\frac{1}{2} \lambda _5 \left[ (\Phi ^\dagger \eta )^2 + \mathrm{h.c.} \right] \nonumber \\&+ \lambda _6 |S|^4 + \lambda _7 |S|^2 |\Phi |^2 + \lambda _8 |S|^2 |\eta |^2. \end{aligned}$$

(2)

After SSB, the scalar fields $\Phi ,~\eta ,~S$ are denoted as

$$\begin{aligned} \Phi&= \left( \begin{array}{cc} G^+ \\ \frac{1}{\sqrt{2}}(v + \varphi + iG^0) \end{array}\right) ,\quad \eta = \left( \begin{array}{cc} \eta ^+ \\ \frac{1}{\sqrt{2}}(\eta _H + i\eta _A) \end{array}\right) ,\quad \nonumber \\ S&= \frac{1}{\sqrt{2}}(v_S + s + iG_S), \end{aligned}$$

(3)

where $v=246$ GeV, $v_S$ is the VEV of S which breaks the $U(1)_{L_\mu -L_\tau }$ symmetry. $G^\pm$, $G^0$ and $G_S$ are corresponding Nambu–Goldstone bosons, which are respectively absorbed by the longitudinal component of the $W^\pm$, Z and $Z^\prime$ gauge bosons. Therefore the VEV $v_S$ provides a mass of $U(1)_{L_\mu -L_\tau }$ gauge boson $Z'$ as $M_{Z'}=g'v_S$. Due to the unbroken $Z_2$ symmetry, the components of inert doublet field $\eta$ ($\eta ^{\pm },~\eta _H,~\eta _A$) do not mix with other scalar fields and their squared masses are simply given as

$$\begin{aligned} M_{\eta ^\pm }^2&= \mu _\eta ^2 +\frac{v^2}{2}\lambda _3 +\frac{v_S^2}{2}\lambda _8, \end{aligned}$$

(4)

$$\begin{aligned} M_{\eta _A}^2&= \mu _\eta ^2 +\frac{v^2}{2}(\lambda _3+\lambda _4-\lambda _5) +\frac{v_S^2}{2}\lambda _8, \end{aligned}$$

(5)

$$\begin{aligned} M_{\eta _H}^2&= \mu _\eta ^2 +\frac{v^2}{2}(\lambda _3+\lambda _4+\lambda _5) +\frac{v_S^2}{2}\lambda _8. \end{aligned}$$

(6)

While for the $Z_2$-even sector, two CP-even scalar components $\varphi$ and s mix with each other. And in the gauge eigenstates $(\varphi ,~s)$, their mass matrix ${\mathcal {M}}_H^2$ is written as

$$\begin{aligned} {\mathcal {M}}_H^2 = \begin{pmatrix} v^2\lambda _1 &{} vv_S\lambda _7 \\ vv_S\lambda _7 &{} 2v_S^2\lambda _6 \end{pmatrix}. \end{aligned}$$

(7)

${\mathcal {M}}_H^2$ can be diagonalized to give the mass eigenstates h and $H_0$. The corresponding squared mass eigenvalues are

$$\begin{aligned} M^2_{h,~H_0}&= \frac{1}{2} \left[ ({\mathcal {M}}_H^2)_{11} + ({\mathcal {M}}_H^2)_{22} \right. \nonumber \\&\quad \left. \pm \sqrt{\left( ({\mathcal {M}}_H^2)_{11} - ({\mathcal {M}}_H^2)_{22}\right) ^2 + 4 \left| {\mathcal {M}}_H^2)_{12}\right| ^2 }\right] , \end{aligned}$$

(8)

The mass and gauge eigenstates are related by

$$\begin{aligned} \left\{ \begin{array}{ll}h &{}= \varphi \cos {\alpha }+s\sin {\alpha },\\ H_0 &{}= -\varphi \sin {\alpha }+s\cos {\alpha },\end{array}\right. \end{aligned}$$

(9)

with the mixing angle given by

$$\begin{aligned} \tan {2\alpha }=\frac{2 ({\mathcal {M}}_H^2)_{12}}{({\mathcal {M}}_H^2)_{11} - ({\mathcal {M}}_H^2)_{22}}. \end{aligned}$$

(10)

Here we assume h as the SM-like Higgs boson with the mass of 125 GeV. Thus, $H_0$ corresponds to an additional singlet-like Higgs boson. In the following, we mainly consider $\alpha \lesssim 0.1$ to avoid various constrains [157]. Finally, the bounded from below condition requires:

$$\begin{aligned}&\lambda _1,~\lambda _2,~\lambda _6 > 0, \end{aligned}$$

(11)

$$\begin{aligned}&\lambda _7 + \frac{1}{\sqrt{2}}\sqrt{\lambda _1\lambda _6}> 0, \quad \lambda _8 + \frac{1}{\sqrt{2}}\sqrt{\lambda _2\lambda _6} > 0, \end{aligned}$$

(12)

$$\begin{aligned}&\lambda _3 + \frac{1}{2}\sqrt{\lambda _1\lambda _6}+\min (0,\lambda _4\pm \lambda _5) > 0. \end{aligned}$$

(13)

At tree level, there is no mixing between Z and $Z'$. But at one-loop level, Z and $Z'$ would mix via the exchange of $\mu ,\nu _\mu ,\tau ,\nu _\tau$ with the loop factor given by [158, 159]

$$\begin{aligned} \Pi ^{\mu \nu }(q^2)= & {} -(q^2g^{\mu \nu }-q^\mu q^\nu )\nonumber \\&\times \left( \frac{g'}{2\pi ^2}\frac{C_V g}{2\cos \theta _W}\right) \int _0^1dx~x(1-x) \text {ln}\nonumber \\&\times \left[ \frac{M_\tau ^2-x(1-x)q^2}{M_\mu ^2-x(1-x)q^2}\right] , \end{aligned}$$

(14)

where $C_V=-\frac{1}{2}+2\sin ^2\theta _W$ and $\theta _W$ is the Weinberg angle. The resulting mixing angle between Z and $Z'$ thus is

$$\begin{aligned} \tan 2\theta _Z=\frac{2\Pi ^{\mu \nu }g_{\mu \nu }}{M^2_{Z'}-M^2_Z}. \end{aligned}$$

(15)

To satisfy the precise measurement of SM Z-boson mass [160], one needs $\tan \theta _Z<10^{-2}$ [159]. For heavy $Z'$ above $M_\tau$, the mixing is suppressed by $M_\tau ^2/M_{Z'}^2$, thus it is easy to satisfy the experimental limit.

2.2 Neutrino mass

The relevant mass terms and Yukawa interactions are flavor dependent

$$\begin{aligned} -{\mathcal {L}}_Y= & {} \frac{1}{2}M_{ee} \overline{N^c_e} N_e +\frac{1}{2}M_{\mu \tau }( \overline{N^c_\mu } N_\tau + \overline{N^c_\tau } N_\mu )\nonumber \\&+ h_{e\mu } (\overline{N^c_e} N_\mu + \overline{N^c_\mu } N_e) S^* + h_{e\tau } (\overline{N^c_e} N_\tau + \overline{N^c_\tau } N_e) S \nonumber \\&+ f_e \overline{L_L^e} {\tilde{\eta }} N_e + f_\mu \overline{L_L^\mu } {\tilde{\eta }} N_\mu + f_\tau \overline{L_L^\tau } {\tilde{\eta }} N_\tau \nonumber \\&+y_e \overline{L_L^e} \Phi e_R + y_\mu \overline{L_L^\mu } \Phi \mu _R + y_\tau \overline{L_L^\tau } \Phi \tau _R+\mathrm{h.c.},\nonumber \\ \end{aligned}$$

(16)

where ${\tilde{\eta }}=(i\sigma _2) \eta ^*$. Hence, after $\Phi (S)$ develops VEV $v(v_S)$, mass matrix of charged lepton $\ell$ and right-handed neutrino $N_\ell$ can be written as

$$\begin{aligned} {\mathcal {M}}_\ell&=\frac{v}{\sqrt{2}}\text {diag}(y_e,y_\mu ,y_\tau ),\quad \nonumber \\ {\mathcal {M}}_N&=\left( \begin{array}{c@{\quad }c@{\quad }c} M_{ee} &{} \frac{v_S}{\sqrt{2}} h_{e\mu } &{} \frac{v_S}{\sqrt{2}} h_{e\tau }\\ \frac{v_S}{\sqrt{2}} h_{e\mu } &{} 0 &{} M_{\mu \tau }e^{i\theta _R} \\ \frac{v_S}{\sqrt{2}} h_{e\tau } &{} M_{\mu \tau }e^{i\theta _R} &{} 0 \end{array} \right) \end{aligned}$$

(17)

With appropriate phase redefinition, all the parameters can be made real, and $\theta _R$ is the CP-violating phase. The symmetric matrix ${\mathcal {M}}_N$ can be diagonalized by an unitary matrix V as

$$\begin{aligned} V^T {\mathcal {M}}_N V= \text {diag}(M_1,M_2,M_3), \end{aligned}$$

(18)

where the lightest one is regard as DM candidate and denoted as N for simplicity in the following discussion.

The neutrino mass is generated at one-loop level via exchanging $\eta$ and $N_\ell$ in the loop. Provided the inert doublet scalar much heavier than the right handed neutrinos, the resulting neutrino mass matrix is approximately given by

$$\begin{aligned} ({\mathcal {M}}_\nu )_{ij}\simeq \frac{\lambda _5 v^2}{8\pi M_0^2} f_i ({\mathcal {M}}_N)_{ij} f_j, \end{aligned}$$

(19)

where $M_0^2=(M^2_{\eta _H}+M^2_{\eta _A})/2$. Similar as ${\mathcal {M}}_N$, the structure of ${\mathcal {M}}_\nu$ corresponds to “Pattern C” of two-zero texture in Ref. [161]. Therefore, only the inverted neutrino mass hierarchy can fit the neutrino oscillation data [65]. Due to the two-zero texture, the nine neutrino parameters are determined by five input parameters. Briefly speaking, the heavy right-handed neutrino mass matrix is determined by $M_{ee}$, $v_S$, $h_{e\mu }$, $h_{e\tau }$ and $M_{\mu \tau }(\theta _R)$. Then for a given ${\mathcal {M}}_N$, the neutrino mass and mixing parameters can be acquired by tuning free parameters $\lambda _5$, $M_0$, $f_e$, $f_\mu$ and $f_\tau$, as long as the following condition is satisfied [66]

$$\begin{aligned} R=\frac{({\mathcal {M}}_\nu )_{12}({\mathcal {M}}_\nu )_{13}}{({\mathcal {M}}_\nu )_{11}({\mathcal {M}}_\nu )_{23}} = \frac{h_{e\mu } h_{e\tau }v_S^2e^{-i\theta _R}}{M_{ee} M_{\mu \tau }}\simeq 0.46\times e^{3.1i}.\nonumber \\ \end{aligned}$$

(20)

3 $R_{K^{(*)}}$ anomaly

In order to account for the $R_K$ and $R_{K^*}$ anomalies, a flavor changing coupling $Z^\prime bs$ is necessary, which is however absent in the original model due to the fact that quarks do not carry $U(1)_{\mu -\tau }$ charges. As a complement, we follow Ref. [162] to extend the original model by introducing a set of heavy vector-like quarks $Q_L \equiv (U_L,~D_L)$, $U_R^c$, $D_R^c$ and their chiral partners ${\tilde{Q}}_{R}\equiv ({\tilde{U}}_R,~{\tilde{D}}_R)$, ${\tilde{U}}_{L}^c$, ${\tilde{D}}_{L}^c$, and whose charge assignment under $SU(2)_L \times U(1)_Y \times U(1)_{\mu -\tau } \times Z_2$ are listed in Table 2.

Table 2 The charge assignment of heavy vector-like quarks under $SU(2)_L \times U(1)_Y \times U(1)_{\mu -\tau } \times Z_2$

The relevant mass terms and the Yukawa interactions of the heavy vector-quarks are

$$\begin{aligned} {\mathcal {L}}_{\text {VLQ}}= & {} S \left( \overline{{{\tilde{D}}}}_R Y_{Qj} P_L d_j + \overline{{{\tilde{D}}}}_L Y_{Dj} P_R d_j \nonumber \right. \\&\left. + \overline{{{\tilde{U}}}}_R Y_{Qj} P_L u_j + \overline{{{\tilde{U}}}}_L Y_{Uj} P_R u_j\right) \nonumber \\&+M_Q \overline{Q}_L {\tilde{Q}}_R + M_D \overline{{{\tilde{D}}}}_L D_R + M_U \overline{{{\tilde{U}}}}_L U_R ~+~ \mathrm{h.c.}~,\nonumber \\ \end{aligned}$$

(21)

After integrating out the heavy vector-like quarks, Eq. (21) induces an effective coupling of $Z^\prime \bar{d}_i d_j$ with the following form

$$\begin{aligned} g^\prime (\varvec{L}^d_{ij} \bar{d}_i \gamma ^\mu P_L d_j Z^\prime _\mu + \varvec{R}^d_{ij} \bar{d}_i \gamma ^\mu P_R d_j Z^\prime _\mu ), \end{aligned}$$

(22)

where the Yukawa-like matrices $\varvec{L}^d_{ij}$ and $\varvec{R}^d_{ij}$ depend on $m_{Q,D}$ and $Y_{Q,D}$ as [162]

$$\begin{aligned} \varvec{L}^d_{ij}= \frac{v^2_S}{2 M_Q^2} (Y_{Qi} Y^*_{Qj}),\quad \varvec{R}^d_{ij} = -\frac{v^2_S}{2 M_D^2}\left( Y_{Di} Y^*_{Dj}\right) . \end{aligned}$$

(23)

Obviously, they are hermitian matrices. In general, these couplings should be complex. For simplicity, all components are taken to be real. We further simplify these matrices by only keeping the flavor-diagonal components and the components which are related to bs transition. Thus, $\varvec{L}^d_{ij}$ and $\varvec{R}^d_{ij}$ take the form

$$\begin{aligned} \varvec{L}^d_{ij} = \frac{v^2_S Y^2_{Q}}{2 M^2_{Q}} \begin{pmatrix} 1 &{}\quad 0 &{}\quad 0\\ 0&{}\quad 1 &{}\quad 1 \\ 0&{}\quad 1 &{}\quad 1 \end{pmatrix},\quad \varvec{R}^d_{ij} = -\frac{v^2_S Y^2_{D}}{2 M^2_{D}} \begin{pmatrix} 1 &{}\quad 0 &{}\quad 0\\ 0&{}\quad 1 &{}\quad 1 \\ 0&{}\quad 1 &{}\quad 1 \end{pmatrix}, \end{aligned}$$

(24)

Considering the best-fit values of Wilson coefficients for $R_{K^{(*)}}$ anomaly in Refs. [163, 164], we found that current data implies $C_9^{\prime \,\mu }\approx 0$. For simplicity, we neglect $C_9^{\prime \,\mu }$ contribution in the later discussions by setting $m_D$ decoupled from the spectrum. After above manipulations, from Eqs. (22) and (24), one obtains following effective Hamiltonian for $b\rightarrow s\mu \mu$ decays

$$\begin{aligned} {\mathcal {H}}^{bs\mu \mu }_{\mathrm{eff}} = \frac{g^{\prime 2} \varvec{L}^d_{23}}{M^2_{Z^\prime }} \bar{s} \gamma ^\mu P_L b {\bar{\mu }}\gamma _\mu \mu ~, \end{aligned}$$

(25)

for heavy $Z^\prime$, and corresponding Wilson coefficient $C_9^\mu$ with muons reads

$$\begin{aligned} C_9^\mu = \frac{g^{\prime ~2}\varvec{L}^d_{23}}{M^2_{Z^\prime }}\frac{\sqrt{2}\pi }{G_F V_{tb}V_{ts}^*\alpha _{em}} = \frac{\pi }{\sqrt{2}G_F V_{tb}V_{ts}^*\alpha _{em}}\left( \frac{Y_{Q}}{M_Q}\right) ^2. \end{aligned}$$

(26)

Taking typical values $G_F = 1.166\times 10^{-5}~\mathrm{GeV}^{-2}$, $V_{tb}V_{ts}^*\approx -4.058\times 10^{-3}$, Eq. (26) yields

$$\begin{aligned}&C_9^\mu \approx -6.43\times 10^9 {\mathrm{GeV}^2} \left( \frac{Y_{Q}}{M_Q}\right) ^2 \nonumber \\&\quad = -1.61\left( \frac{Y_{Q}^2 }{0.1}\right) \left( \frac{20 \mathrm{TeV}}{M_Q}\right) ^2~. \end{aligned}$$

(27)

To interpret the $R_{K^{(*)}}$ anomaly, the $1\sigma$ range $C_9^\mu \in [-1.10,-0.79]$ is required [164]. It is clear that the coefficient $C_9^\mu$ only depends on parameters $Y_Q$ and $M_Q$, leaving the $U(1)_{L_\mu -L_\tau }$ parameters $m_{Z'}$ and $g'$ free to choose. In the following discussion, we fix $Y_Q=0.1215$ and $M_Q=10~\mathrm{TeV}$ to acquire the best fit value of $C_9^\mu \approx -0.95$ [164].

3.1 Constraints

Although interpretation of $R_{K^{(*)}}$ in above section do not depends on $M_{Z'}$ and $g'$, the relevant parameter space of $Z'$ is constrained by the following experiments:

• Moun $g-2$ and neutrino trident production In our model, the $Z^\prime$ contribution to the muon magnetic moment anomaly is given as

  $$\begin{aligned} \Delta {a_\mu }&\equiv \frac{(g-2)_\mu }{2}\nonumber \\&=\frac{g^{\prime 2}}{8\pi ^2}\int ^1_0 \frac{2x^2 (1-x)}{x^2+(M_{Z^\prime }/M_\mu )^2(1-x)}dx, \end{aligned}$$

  (28)

  However, the allowed parameter space is tightly constrained by neutrino trident production [165], i.e., $\nu _\mu N\rightarrow \nu _\mu N\mu ^+\mu ^-$ process. In the heavy $Z^\prime$ case, the normalized cross section expressed as [162]

  $$\begin{aligned} \frac{\sigma }{\sigma _{\mathrm{SM}}} \simeq \frac{1+\big (1+4s_W^2 + 8 \frac{g^{\prime 2}}{g^2_2}\frac{M^2_W}{M^2_{Z^{\prime }}}\big )^2}{ 1+(1+4s_W^2 )^2}. \end{aligned}$$

  (29)

  In this paper, we consider the CCFR measurement $\sigma /\sigma _{\mathrm{SM}}=0.82\pm 0.28$ [166].

• $B_s-\bar{B}_s$ mixing In our model, the flavor changing vortex $Z^\prime s b$ leads to tree level $B_s$ mixing via $Z^\prime$ exchange. In addition, Yukawa interactions in Eq. (21) also contribute to $B_s$ mixing via box diagram, where singlet scalar and heavy vector-like quark running in the loop. The modification of mixing amplitude for heavy $Z^\prime$ yield [162]

  $$\begin{aligned} \frac{\Delta {M_{12}}}{M_{12}^\text {SM}}\simeq & {} \left( \frac{Y_Q}{M_Q}\right) ^4 \left[ \left( \frac{M_{Z^\prime }}{g^\prime }\right) ^2+\frac{M^2_Q}{16\pi ^2}\right] \nonumber \\&\times \left[ \frac{g^4_2}{16\pi ^2}\frac{1}{M^2_W}(V_{tb}V^*_{ts})^2S_0\right] ^{-1}, \end{aligned}$$

  (30)

  with $S_0\approx 2.3$. Current experimental measurement allows $\Delta M_{12}/ M_{12}^{\text {SM}}\lesssim 15\%$ [167].

• Branching ratio for $t\rightarrow cZ^\prime$ This is also induced by the left-handed $t\rightarrow c$ current, which is related to the left-handed $b\rightarrow s$ current by SU(2)$_L$ symmetry. Provided $Y_{Qt}\sim Y_{Qb}$, $Y_{Qc}\sim Y_{Qs}$ [168], the branching ratio is

  $$\begin{aligned} \mathrm{BR}(t\rightarrow cZ^\prime )_{\mathrm{LH}}&\simeq \frac{(1-x_{Z^\prime })^2(1+2x_{Z^\prime })v^2}{8(1-x_W)^2(1+2x_W)}\nonumber \\&\quad \left( \frac{Y_Q}{M_Q}\right) ^4 \left( \frac{M_{Z^\prime }}{g^\prime }\right) ^2, \end{aligned}$$

  (31)

  where $x_W\equiv M_W^2/M_t^2$, $x_{Z^\prime } \equiv M_{Z^\prime }^2/M_t^2$ and we have set $M_c^2/M_t^2$, $M_b^2/M_t^2 \rightarrow 0$. The decay $t\rightarrow cZ^\prime$ followed by $Z^\prime \rightarrow \ell ^+\ell ^-$ ($\ell =\mu ,\tau$) can be searched for in $t{\bar{t}}$ events at the LHC. It is similar to $t\rightarrow qZ$ ($q=u,c$) decay, which has been searched for by the ATLAS [169] and CMS [170] experiments using $t{\bar{t}} \rightarrow Zq+Wb$ with leptonically decaying Z and W, resulting in a final state with three charged leptons. Reinterpreting the CMS limits for $t\rightarrow cZ$ to the case for $t\rightarrow cZ^\prime$ by a simple scaling of Z and $Z^\prime$ decay branching ratios into the charged leptons ($\ell =e,\mu$), Ref. [168] found $\mathrm{BR}(t\rightarrow cZ')\lesssim 10^{-4}$.

• Branching ratio for $Z\rightarrow 4\ell$ In our model, the $Z\rightarrow 4\ell$ decay will receive a significant contribution from $Z\rightarrow \mu ^+\mu ^-Z^\prime$ followed by $Z^\prime \rightarrow \mu ^+\mu ^-$ for the case of $M_{Z^\prime }< M_{Z}$. The ATLAS [171] and CMS [172, 173] collaborations both have set upper limits on the branching fraction of the Z boson decay to four charged leptons. In particular ATLAS has set an upper limit on $\mathrm{BR}(Z\rightarrow 4\mu )=(4.2\pm 0.4)\times 10^{-6}$ with the combined 7 TeV and 8 TeV dataset [171]. Using $77.6~\mathrm{fb}^{-1}$ data at 13 TeV, CMS recently sets a more stringent upper limits of $10^{-8}\sim 10^{-7}$ on the branching ratio BR$(Z\rightarrow Z'\mu \mu )$BR$(Z'\rightarrow \mu \mu )$ [172]. In this work, we adopt the dedicated limits on $L_\mu -L_\tau$ model provided by CMS [172].

• LHC $Z^\prime$ constraints on dilepton final state. In our model, $Z^\prime$ boson will be produced at LHC through the flavor conserving process $q\bar{q}\rightarrow Z^\prime$ and flavor violating process $b\bar{s}\rightarrow Z^\prime$ (and its conjugate process). Therefore, searches of heavy resonance in the dimuon final state by ATLAS [174, 175] and CMS [176] tightly constrain the parameter space. In particular, ATLAS [174, 175] has set a $95\%$ C.L. upper limit on $\sigma (pp\rightarrow Z^\prime + X) \mathrm{BR}(Z^\prime \rightarrow \mu ^+\mu ^-)$ in the 150 GeV $\lesssim M_{Z^\prime } \lesssim$ 5 TeV mass range, with the 13 TeV and $\sim$ 13 fb$^{-1}$ dataset.

Other experiments, such as $\tau$ decays, are less strict than the above ones [162], we thus do not take into account in this paper. The above mentioned constraints are shown in Fig. 3.

Fig. 1

Distribution of dominant annihilation channels in the $g'$-$M_Z'$ plane. Survived samples in left panel satisfy constraints from relic density, while those in right panel satisfy constraints from both relic density and direct detection

4 Dark matter phenomenology

In this section, we further investigate the phenomenology of Majorana fermion DM N for high mass $Z'$. The model is implemented in FeynRules [177]. The calculation of DM relic density and DM-nucleon scattering cross section are performed with the help of micrOMEGAs [178]. The possible annihilation channels in this gauged $U(1)_{L_\mu -L_\tau }$ scotogenic model are listed in the following.

• $NN\rightarrow \ell ^+\ell ^-, \nu _\ell \nu _\ell (\ell =e,\mu ,\tau )$ is mediated by the inert scalar doublet $\eta$ via the Yukawa coupling $f_\ell$. Tightly constrained by LFV, $f_\ell \lesssim 0.01$ is usually needed for electroweak scale N and $\eta$ [40,41,42], thus contributions of this channel are negligible.

• $NN\rightarrow Z^{'*}\rightarrow \ell ^+\ell ^-, \nu _\ell \nu _\ell (\ell =\mu ,\tau )$ via the gauge coupling $g'$ provides a new s-channel process for DM annihilation. Different from the gauged $U(1)_{B-L}$ case [43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59], this channel exclusively generates muon, tau leptons and neutrinos.

• $NN\rightarrow h^*/H_0^*\rightarrow W^+W^-,b\bar{b},\ldots$ via the Yukawa coupling $h_{e\mu }$, $h_{e\tau }$ is also a s-channel process. Previous study neglected this channel by assuming tiny mixing angle $\alpha$ [66]. In this paper, a not too small mixing angle $\alpha$ is considered.

• $NN\rightarrow Z' Z', Z'H_0, H_0H_0, hH_0$ are also possible if kinematically allowed. As shown latter, the $NN\rightarrow Z'Z', Z'H_0(\rightarrow Z'Z')$ channel is possible to interpret the AMS-02 positron excess.

Fig. 2

Same as Fig. 1, but in the $M_{Z'}$-$M_N$ plane

To illustrate the effects of above various annihilation processes, we implement a random scan over the following parameter space

$$\begin{aligned} g'\in [0.001,1]&,&M_{Z'}\in [10,5000]~\mathrm{GeV}, \nonumber \\ \alpha \in [0.01,0.1]&,&M_{H_0} \in [0,\sqrt{4\pi } M_{Z'}/g'], \nonumber \\ h_{e\mu ,e\tau }\in [0,4\pi ]&,&M_{ee,\mu \tau }\in [10,5000]~\mathrm{GeV}, \end{aligned}$$

(32)

and assign the dominant annihilation channel to the survived samples under constraints from relic density and direct detection. For relic density, we use the combined Planck+WP+highL+BAO $2\sigma$ value, i.e., $0.1153<\Omega h^2<0.1221$ [179]. As for direct detection, we adopt the combined limits provided by XENON1T [129, 130] and PandaX-II [131]. Meanwhile, to satisfy the observed neutrino oscillation parameters, we further require |R| defined in Eq. (20) in the interval [0.4, 0.5] with $\theta _R=\pi$ for simplicity [66].

Due to the Majorana nature of N DM, the DM-nucleon scattering cross section mediated by $Z'$ is suppressed, and is actually dominant by Higgs exchange. In this way, the spin-independent cross section is given by [180]

$$\begin{aligned} \sigma ^\text {SI}=\frac{h_N^2 \mu _N^2 M_n^2 f_n^2}{2 \pi v^2} \sin ^2 2\alpha \left( \frac{1}{M_h^2}-\frac{1}{M_{H_0}^2}\right) ^2, \end{aligned}$$

(33)

where $h_N=(V_{eN}V_{\mu N}+V_{\mu N}V_{e N}) h_{e\mu }+ (V_{eN}V_{\tau N}+V_{\tau N}V_{e N}) h_{e\tau }$ is the effective DM-S coupling, $M_n\approx 0.939~\mathrm{GeV}$ is the averaged nucleon mass, $\mu _N=M_n M_N/(M_n+M_N)$ is the DM-nucleon reduced mass, $f_n\approx 0.345$ is the nucleon matrix element. Clearly, the cross section is proportional to $\sin ^22\alpha$, therefore a smaller mixing angle $\alpha$ is also preferred by direct detection.

In Figs. 1 and 2, we depict the distribution of survived samples in the $g'$-$M_{Z'}$ and $M_{Z'}$-$M_N$ plane respectively. From Fig. 1, we aware that correct relic density could be realized with $g'\gtrsim 0.02$ and $M_{Z'}\gtrsim 20~\mathrm{GeV}$ during our scan, but the direct detection would exclude those points with $M_{Z'}\lesssim 100~\mathrm{GeV}$. Note that a little points dominant by $hH_0$ is survived under relic density, but such points are fully excluded by direct detection. Distributions of different annihilation channels in the $M_{Z'}$-$M_N$ plane as shown in Fig. 2 are clearer for different annihilation channels. The $NN\rightarrow Z'^*\rightarrow \ell ^+ \ell ^-$ channel is dominant in the region $M_{Z'}\sim 2 M_{N}$. When $M_{Z'}<M_N$ or $M_{H_0}<M_N$, the dominant channels become $NN\rightarrow Z'Z'$, $ZH_0$ and $H_0H_0$. For the s-channel Higgs portal dominant, the $NN\rightarrow b\bar{b}$ channel needs $M_N\sim M_h/2\approx 60~\mathrm{GeV}$, while $NN\rightarrow W^+W^-$ channel requires $M_N\gtrsim 500~\mathrm{GeV}$.

Fig. 3

Survived samples in the $g'$-$M_Z'$ plane with constraints from Sect. 3.1. The constraint $t\rightarrow c Z'$ excludes $g'\lesssim 4\times 10^{-4}$, which is too small to show here. The cyan circle points satisfy relic density only, while the orange triangle points further satisfy direct detection. Blue circle points satisfy relic density and neutrino oscillation, and the red triangle points further satisfy direct detection

In Fig. 3, we show combined results from relic abundance, direct detection, neutrino oscillation as well as various constraints on $Z'$ discussed in Sect. 3.1. The neutrino trident production process has excluded the $(g-2)_\mu$ favor region, and set the most stringent upper limit on $g'$ for $M_{Z'}\gtrsim 50~\mathrm{GeV}$. For $M_{Z'}\lesssim 50~\mathrm{GeV}$, the most stringent upper limits comes from $Z\rightarrow 4\mu$ search at LHC. Meanwhile, the $B_s$ mixing has set an lower limit on $g'$. It is clear that a few survived red triangle samples are not excluded by neutrino trident production and $B_s$ mixing. Hence, viable parameter space is obtained to explain DM, neutrino mass and $R_{K^{(*)}}$ anomaly simultaneously.

5 AMS-02 positron excess

5.1 AMS-02 positron flux

In this section, we discuss the AMS-02 positron excess and relevant constraints from indirect detections. Recently, the AMS Collaboration has released latest result of positron spectrum that extend the maximal measurement energy up to 1 TeV [181], which is used in our fitting. For a given model parameters in Eq. (32), the positron flux can be expressed as

$$\begin{aligned} \Phi _{e^+}(f_{e^+},\phi _{e^+}^\odot , BF) = f_{e^+} \Phi _{e^+}^{\mathrm{bkg}}(\phi _{e^+}^\odot )+\Phi _{e^+}^{\mathrm{DM}}(\phi _{e^+}^\odot , BF), \end{aligned}$$

(34)

where $f_{e^+}$ is the normalization factors which take into account the uncertainty of the astrophysical background and varying in the range [0,  5] in our fitting. The fluxes of charged CR particles are periodically modulated according to the solar activity due to their interactions with the heliosphere magnetic field. The modulation is more important for low energy CR particles and can be described by using the force field approximation [182]. In this approximation, the modulated spectrum $\Phi _{\mathrm{mod}}(E_k)$ and unmodulated one $\Phi (E_k)$ is related by following formula:

$$\begin{aligned} \Phi _{\mathrm{mod}}(E_k)=\frac{(E_k+M)^2-m_{\mathrm{CR}}^2}{(E_k+m_{\mathrm{CR}}+|e|\phi ^\odot )^2-m_{\mathrm{CR}}^2} \Phi (E_k+|e|\phi ^\odot ), \end{aligned}$$

(35)

where $m_{\mathrm{CR}}=m_e$ ($m_p$) for position (antiproton) flux, $\phi ^\odot$ is the modulation potential, and $E_k$ is the observed kinetic energy. Note that the force field approximation is an over simplified model. The modulation potential $\phi ^\odot$ is just an effective parameter that indicate the total effect in the solar modulation. In fact, different CR particles would always require different modulation potentials. We therefore consider the modulation potentials of positron and antiproton as two free parameters in the fitting.

We calculate the positron and antiproton flux resulted from DM annihilation using micrOMEGAs [178]. For the DM density distribution in the galactic halo, we have used the Navarro–Frenk–White (NFW) density profile [183] with the local density $\rho _\odot =0.4\,\mathrm {\,GeV\,cm}^{-3}$. The background fluxes are obtained by solving the diffusion equation for cosmic-ray particles using a widely used galactic cosmic-ray propagation model. To take into account the convection/reacceleration effect and the complex electron energy losses during the diffusion, we adopt the public code GALPROPv54 [184, 185]. The relevant parameters for these process, such as the diffusion coefficient and the convection velocity, ought to be determined by fitting to the B/C and proton data. Here we choose the these parameters following the diffusion + convection (DC) case in Ref. [186] to derive both the secondary positron and antiproton fluxes. In all, we set the diffusion coefficient $D(R)=1.95\times 10^{28}(R/4.71\mathrm {\,GV})^{0.51}\mathrm {\,cm^2\,s^{-1}}$, the gradient of covection velocity $\mathrm {d}V/\mathrm {d}z=4.2\mathrm {\,km\,s^{-1}\,kpc^{-1}}$ and the proton injection with a power-index 2.336. The $\chi ^2$ function is defined as

$$\begin{aligned}&\chi ^2_{e^+}(f_{e^+},\Phi _{e^+}^\odot , BF)\nonumber \\&\quad =\sum _i\frac{\left[ \Phi _{e^+,i}(f_{e^+},\Phi _{e^+}^\odot , BF)-\Phi _{e^+,i}^{\mathrm{AMS}}\right] ^2}{(\sigma _{e^+,i}^{\mathrm{AMS}})^2}, \end{aligned}$$

(36)

where i runs over all the data points. $\Phi _{e^+,i}$ and $\sigma _{e^+,i}^{\mathrm{AMS}}$ are respectively the relevant observables (positron fraction in this case) and corresponding experimental errors (stat+syst) taken from [136, 137].

It has been known that in order to fit AMS-02 positron data, a large enhancement with $BF\sim \mathcal{O}(10^3)$ is required for annihilation cross section in the Galaxy with $v\sim 10^{-3}$ than that in the freeze-out temperature with $v\sim 10^{-1}$. In addition, the annihilation final states should be leptophilic to avoid antiproton constraint [186]. Even so, this scenario is still challenged by limits from extragalactic $\gamma$-ray background [187] and CMB observations [179, 188,189,190]. We will discuss these constraints in more detail in section 5.2. Here we first illustrate how to obtain a large BF in a consistent way in our model. The two common methods widely used to simultaneously realize the correct relic abundance and a large BF are so-called Breit–Wigner mechanism [191,192,193] and Sommerfeld enhancement [194,195,196]. In the former case, two DM particles annihilate via the s-channel exchange of a heavy mediator, then the annihilation cross section are resonantly enhanced when mediator mass is close to twice of DM mass. Notably, it has been shown that Breit–Wigner mechanism can potentially relax the tension between positron excess and CMB observations due to the evolution of velocity dependent annihilation cross section at different cosmic epochs (freeze-out, recombination and present) [197]. Unfortunately, this mechanism has less effect on our model. Based on the discussion in Sect. 4, the only important s-channel annihilation is $NN\rightarrow Z^{'*}\rightarrow \ell ^+\ell ^-$, which is p-wave suppression in the Galaxy since N is Majorana DM. As a consequence, annihilation cross section is not large enough even in the resonance regions.

Fig. 4

The benchmarks in our model to fit AMS-02 positron data. Here cyan, orange and purple diamonds (red ,green and blue stars) respectively corresponding to $M_N =1,~1.5,~2$ TeV for $NN\rightarrow Z^\prime Z^\prime$ ($NN\rightarrow Z^\prime H_0$) annihilation channel. For completeness, various limits in Fig. 3 are shown

We therefore focus on Sommerfeld enhancement in the s-wave DM annihilation. This mechanism is due to the loop correction of annihilation cross section with exchange infinite number of vector or scalar mediators. In the non-relativistic limit, the velocity-dependent correction to DM annihilation can be computed by numerically solving the radial Schrödinger equation with the attractive spherically symmetric Yukawa potential $V(r)=-\alpha ^\prime e^{-M^\prime }/r$, here $\alpha ^\prime$ and $M^\prime$ respectively denote coupling and mass of mediator. The Sommerfeld enhancement factor $S_E$ then evaluated by the radial wave function at the origin, $|\psi _k(0)|^2$. In following, we will use the semi-analytic formula introduced in Refs. [198, 199] for a illustration. Both $Z^\prime$ and $H_0$ can serve as mediators in our model, with corresponding annihilation channel are respectively $NN\rightarrow Z^\prime Z^\prime \rightarrow 2\ell ^+2\ell ^-+2\nu _\ell 2\nu _\ell$ and $NN\rightarrow Z^\prime H_0(\rightarrow Z'Z') \rightarrow 3\ell ^+3\ell ^-+3\nu _\ell 3\nu _\ell$ with $\ell =\mu ,\tau$. However, annihilation channel $NN\rightarrow Z^\prime Z^\prime$ is difficult to give desired BF through Sommerfeld enhancement and leaving $NN\rightarrow Z' H_0$ as only available channel, which is due to the fact that for $NN\rightarrow Z' Z'$ channel, its Sommerfeld enhancement factor determined by relic abundance and by positron excess are incompatible.

Table 3 The DM information for benchmarks of $NN\rightarrow Z^\prime Z^\prime$ annihilation channel. Here we choose $h_N=0.5$ and $\langle \sigma v\rangle _0$ (in units of $\mathrm{cm}^3~\mathrm{s}^{-1}$) denotes the thermally averaged DM annihilation cross section at freeze out, $\langle \sigma v\rangle _{BF}\equiv BF\times \langle \sigma v\rangle _0$ denotes annihilation cross section in the Galactic halo which required by interpreting AMS-02 positron excess, and $\langle \sigma v\rangle _{\mathrm{CMB}}$ the annihilation cross section limited by CMB observation. All of masses in units of $\mathrm{GeV}$ and solar modulation in units of MV

Fig. 5

The Sommerfeld enhancement factor $S_E$ in our model as a function of mediator mass $M_{Z^\prime }$ (left) and $M_{H_0}$ (right) for the benchmark parameters in Tables 3 and 4

To see this, we plot contours of correct relic abundance on the $g'-M_{Z'}$ plane for benchmark masses $M_N=1,~1.5,~2$ TeV in Fig. 4 with combine various constraints from Fig. 3. Since annihilation cross section of $NN\rightarrow Z' Z'$ channel scales as $g'^4M^2_N/M^4_{Z'}$, relic abundance is entirely fixed by ratio $g'/M_{Z'}$ for a given $M_N$. Thus three curves presented just as straight line in the figure, and $Z'$ masses below 50 GeV are excluded by $Z\rightarrow 4\mu$ LHC direct search. We then select largest $(g',~M_{Z'})$ point in each curve as benchmark (cyan, orange and purple diamonds), their values are listed in Table 3. Corresponding $S_E$ curves for three benchmark $g'$ are shown in left panel of Fig. 5. Comparing the fitting values of $\langle \sigma v\rangle _{BF}\equiv BF\times \langle \sigma v\rangle _0$ in table 3 with $S_E$ values in the figure for the same $M_{Z'}$, it obviously fails to satisfy the requirement $S_E\simeq BF$. We further evaluate all of $(g',~M_{Z'})$ in relic abundance curves, none of them can match above condition. It boils down to the fact that both relic abundance and Sommerfeld enhancement factor $S_E$ share the same coupling $g'$ for given $M_N$ and $M_{Z'}$, which is hardly to tune its value to satisfy two requirements simultaneously.

As a consequence, we appeal to another annihilation channel $NN\rightarrow Z^\prime H_0$. This channel has advantage that its relic abundance and $S_E$ depend on parameters $(g',~h_N,~M_{Z'}, M_{H_0})$, which does not suffer from tight correlation between relic abundance and $S_E$ with such more degree of freedom. Similarly, we presented three benchmarks of $NN\rightarrow Z^\prime H_0$ channel in Fig. 3 (red, green and blue stars) and in Table 4 for the same $M_N$. Here we take $g^\prime =3\times 10^{-3}$ which is compatible with the LHC $Z\rightarrow 4\mu$ direct search limit for light $Z^\prime$, and $M_{H_0}$ for each benchmark has been chosen such that $S_E\simeq BF$, as is shown in the right panel of Fig. 5.

The positron flux predicted by benchmarks in Tables 3 and 4 are shown in Fig. 6 with AMS-02 data. We found that both $NN\rightarrow Z' Z'$ and $NN\rightarrow Z^\prime H_0$ annihilation channels can provide good fitting for $M_N$ in the range of 1–1.5 TeV, while 2 TeV benchmark results in too much excess. Despite only $NN\rightarrow Z^\prime H_0$ channel gives consistent interpretation in Sommerfeld enhancement scenario, from the phenomenological viewpoint, we still treat $NN\rightarrow Z' Z'$ channel as a valid candidate. The relevant constraints for both two channels are investigated in the next section.

Table 4 Similar with Table 3, but for $NN\rightarrow Z^\prime H_0$ annihilation channel. Here we set $g^\prime =3\times 10^{-3}$ which is compatible with the LHC $Z\rightarrow 4\mu$ direct search limit for light $Z^\prime$ boson. The values of $M_{H_0}$ have been chosen such that the resulted Sommerfeld enhancement factors are match to fitted boost factors, i.e., $S_E\simeq BF$

Fig. 6

The positron fluxes predicted by benchmarks of $NN\rightarrow Z'Z'$ (Table 3) and $NN\rightarrow Z'H_0$ (Table 4) annihilation channels with the fitting results of AMS-02 data

5.2 Constraints from other indirect detections

5.2.1 AMS-02 antiproton constraint

Given the benchmarks to explain the positron excess, we now exam the constraints from other indirect detections. As we mentioned in Sect. 5.1, the relevant limits come from antiproton flux, EGRB measurement and impact of energy deposition on CMB anisotropy. From Eq. (22), the effective coupling of $Z^\prime \bar{d}_i d_j$ leads to antiproton flux. Although highly suppressed by Yukawa-like matrices $\varvec{L}^d_{ij}$ and $\varvec{R}^d_{ij}$, we still need to investigate whether the predict antiproton flux conflicts with current observation. Similar with Eqs. (34) and (36), the antiproton flux and $\chi ^2$ function are respectively given by

$$\begin{aligned} \Phi _{\bar{p}}(f_{\bar{p}},\phi _{\bar{p}}^\odot , BF)= & {} f_{\bar{p}} \Phi _{\bar{p}}^{\mathrm{bkg}}(\phi _{\bar{p}}^\odot )+\Phi _{\bar{p}}^{\mathrm{DM}}(\phi _{\bar{p}}^\odot , BF),\nonumber \\ \chi ^2_{\bar{p}}(f_{\bar{p}},\Phi _{\bar{p}}^\odot , BF)= & {} \sum _i\frac{\left[ \Phi _{\bar{p},i}(f_{\bar{p}},\Phi _{\bar{p}}^\odot , BF)-\Phi _{\bar{p},i}^{\mathrm{AMS}}\right] ^2}{(\sigma _{\bar{p},i}^{\mathrm{AMS}})^2}.\nonumber \\ \end{aligned}$$

(37)

$f_{\bar{p}}$ is normalization factors for antiproton background and also set to vary in the range [0,  5]. Note that $f_{\bar{p}},~\phi _{\bar{p}}^\odot$ should be different to $f_{e^+},~\phi _{e^+}^\odot$ since the astrophysical sources and propagation processes are distinct for positron and antiproton. On the other hand, BF should be the same due to the fact that enhancement of annihilation cross section is universal for lepton and quark final states. The resulted antiproton flux for benchmarks to interpret positron excess are plotted in Fig. 7 with AMS-02 measurement [137], and best fit parameter values and $\chi ^2_{\mathrm{min}}$ are listed in Tables 3 and 4. All of benchmarks just have same values for parameters and $\chi ^2_{\mathrm{min}}$, which means that antiproton flux from DM contribution is so small that $\chi ^2$ is entirely determined by background. We thus conclude that our model is totally safe for antiproton constraint.

Fig. 7

The antiproton flux predicted by benchmarks in Tables 3 and 4 with the AMS-02 data

5.2.2 Fermi-LAT EGRB constraint

The next important constraint we consider comes from the EGRB measured by Fermi-LAT collaboration [187]. For calculation of EGRB flux, we follow the procedures in Ref. [200, 201]. The flux of at redshift z is given as

$$\begin{aligned} \frac{d\Phi _{\mathrm{EGB}}}{dE_\gamma }&= \frac{c(1+z)^2{\bar{\rho }}_{\mathrm{DM}}^2BF\langle \sigma v\rangle _0}{8\pi M^2_N}\nonumber \\&\quad \times \int _z^\infty dz'\frac{(1+z')^3B(z',m_{\mathrm{min}})}{H(z')}\frac{dN}{dE_\gamma '}e^{-\tau (z,z',E_\gamma ')}, \end{aligned}$$

(38)

where $E_\gamma '=E_\gamma (1+z')/(1+z)$, ${\bar{\rho }}_{\mathrm{DM}}=\rho _c\Omega _{\mathrm{DM}}$ is the average density of DM with $\rho _c$ the critical density of the Universe at present, $H(z)=H_0\sqrt{(\Omega _{\mathrm{DM}}+\Omega _b)(1+z)^3+\Omega _\Lambda }$ is the Hubble function.

The $\gamma$-ray generation spectrum for per DM annihilation, $dN/dE_\gamma '$, is dominated by two components:

$$\begin{aligned} \frac{dN}{dE_\gamma '}=\left. \frac{dN}{dE_\gamma '}\right| _{\mathrm{FSR}}+\left. \frac{dN}{dE_\gamma '}\right| _{\mathrm{IC}}, \end{aligned}$$

(39)

where $\left. \frac{dN}{dE_\gamma '}\right| _{\mathrm{FSR}}$ corresponding to $\gamma$-rays produced from the final state radiation (FSR) of primary charged lepton final states ($\mu$ and tau in our model) due to DM annihilation. $\left. \frac{dN}{dE_\gamma '}\right| _{\mathrm{IC}}$ characterizes the $\gamma$-rays resulted from Inverse Compton (IC) scattering between secondary electrons/positrons and CMB photons. For the FSR photon spectrum, we adopt analytical formulas in Refs. [202] and [203]

$$\begin{aligned} \left. \frac{dN}{dx}\right| ^\mu _{\mathrm{FSR}}= & {} \frac{\alpha _{\mathrm{e.m.}}}{\pi }\frac{1+(1-x)^2}{x}\ln \left( s(1-x)/m^2_\mu \right) \,,\nonumber \\ \left. \frac{dN}{dx}\right| ^\tau _{\mathrm{FSR}}= & {} x^{-1.31}(6.94x-4.93x^2-0.51x^3)e^{-4.53x}\,,\nonumber \\ \end{aligned}$$

(40)

where $\alpha _{\mathrm{e.m.}}$ is the fine-structure constant, $s=4M^2_N$ and $x=E/M_N$. The $\gamma$-ray photons from IC component is given as [204]

$$\begin{aligned} \left. \frac{dN}{dE}\right| _{\mathrm{IC}}=\int d\epsilon n_\gamma (\epsilon ) \int dE_e\frac{dn}{dE_e}\sigma _{\mathrm{KN}}(\epsilon ,E_e,E)\,. \end{aligned}$$

(41)

In above equation, $n_\gamma (\epsilon )$ is the photon number density of background radiation, and $\sigma _{\mathrm{KN}}(\epsilon ,E_e,E)$ is the Klein-Nishina cross section. $dn/dE_e$ is the electron/positron energy spectrum after propagating, which is related to production spectrum $dN_e/dE'_e$ by following equation

$$\begin{aligned} \frac{dn}{dE_e}=\frac{1}{b(E_e,z)}\int ^{m_N}_{E_e}dE'_e\frac{dN_e}{dE'_e}\,, \end{aligned}$$

(42)

with $b(E_e,z)\simeq 2.67\times 10^{-17}(1+z)^4)(E_e/\mathrm{GeV})^2\mathrm{GeV}\mathrm{s}^{-1}$ is the energy loss rate [204]. Eq. (38) contains a cosmological boost factor which account for the effect of DM halo clustering, $B(z)\equiv \langle (1+\delta (z))^2\rangle =1+\langle \delta ^2(z)\rangle$. We here adopt a halo model that approximates the matter distribution in the Universe as a superposition of DM halos and $B(z,m_{\mathrm{min}})$ can be expressed as

Fig. 8

Cosmological boost factor B(z) as a function of redshift for Maccio concentration model with $M_{\mathrm{min}}=10^{-6}~M_{\bigodot }$

Fig. 9

Comparison of the EGRB flux produced by benchmarks in Tables 3 and 4 with the Fermi-LAT measurements

$$\begin{aligned} B(z,M_{\mathrm{min}})=1+\frac{\Delta _c}{3{\bar{\rho }}_{m,0}}\int ^\infty _{M_{\mathrm{min}}}dM M \frac{dn}{dM}(M,z)f[c(M,z)]. \end{aligned}$$

(43)

where ${\bar{\rho }}_{m,0}$ is the matter density at present, $\Delta _c\simeq 200$ is the overdensity at which the halos are defined and $M_{\mathrm{min}}$ is the minimal halo mass used in integration. $\frac{dn}{dM}(M,z)$ is the halo mass function with the universal form

$$\begin{aligned} \frac{dn}{dM}(M,z)=\frac{{\bar{\rho }}_{m,0}}{M^2}\nu f(\nu )\frac{d\log \nu }{d\log M}. \end{aligned}$$

(44)

In above equation, the parameter $\nu =[\delta _c(z)/\sigma (M)]^2$. $\delta _c(z)$ is the critical overdensity and the $\sigma (M)$ is the variance of the linear density field in spheres containing a mean mass M. c(M, z) represents the halo concentration parameter function and the function f(c) for the halos with the NFW density profile is given as

$$\begin{aligned} f(c)=\frac{c^3}{3}\left[ 1-\frac{1}{(1+c)^3}\right] \left[ \log (1+c)-\frac{c}{1+c}\right] ^{-2}. \end{aligned}$$

(45)

We choose Maccio concentration model [205] and set $M_{\mathrm{min}}=10^{-6}~M_{\odot }$ in calculation, and evolution of $B(z,m_{\mathrm{min}})$ with redshift is presented in Fig. 8. In last, $\tau (E_\gamma ', z,z',)$ is the the optical depth of $\gamma$-ray photon with energy $E'$ and propagating from $z'$ to z. Which can be expressed as

$$\begin{aligned} \tau (E_\gamma ',z,z')=\int ^{z'}_{z}dz''\frac{\alpha (E_\gamma '',z'')}{H(z'')(1+z'')}, \end{aligned}$$

(46)

where $\alpha (E_\gamma ,z)$ is the absorption coefficient and $E_\gamma ''=E_\gamma '(1+z'')/(1+z')$. For detailed description of interactions and corresponding absorption coefficients for photon propagation which are taken into account, see Refs. [200, 201].

In Fig. 9, we show the total EGRB flux for benchmarks in Tables 3 and 4, with the Fermi-LAT measurement [187]. We found that our benchmarks are marginally compatible with observation. However, for different concentration model and smaller minimal halo, they could be potentially ruled out.

5.2.3 Planck CMB constraint

The last constraint necessarily need to consider is the effect of DM annihilation on CMB anisotropy. Annihilation of DM to SM particles between recombination and reionization epoch can inject and deposit energy into intergalactic medium (IGM) through produced electrons, positrons and photons via photoionization, Coulomb scattering, Compton processes, bremsstrahlung and recombination. The primary effect of these processes are to alter ionization fraction and left an imprint on spectrum of CMB anisotropy. The injection power into the IGM per unit volume at redshift z is given by as [189]

$$\begin{aligned} \left( \frac{dE}{dVdt}\right) _{\mathrm{injected}}=\rho _{\mathrm{DM,0}}^2(1+z)^6\frac{g \langle \sigma v\rangle }{M_{\mathrm{DM}}}, \end{aligned}$$

(47)

where $\rho _{\mathrm{DM,0}}$ is the present DM density, and the degeneracy $g=1$ for our Majorana DM N. the relationship between deposited energy and injected one can be parameterized as

$$\begin{aligned} \left( \frac{dE}{dVdt}\right) _{\mathrm{deposited}}=f_{\mathrm{eff}}(z)\left( \frac{dE}{dVdt}\right) _{\mathrm{injected}}, \end{aligned}$$

(48)

where $f_{\mathrm{eff}}(z)$ denotes dimensionless efficiency factor. It is conventional to define function $p_{\mathrm{ann}}(z)=f_{\mathrm{eff}}(z)\langle \sigma v\rangle /M_{\mathrm{DM}}$, which contains full information about the CMB constraint. In specific, for a given primary annihilation final state i with the annihilation spectrum of positron $dN^i_{e^+}/dE$ and photon $dN^i_\gamma /dE$, $f_{\mathrm{eff}}(z)$ can be weighted as [190],

$$\begin{aligned}&f^i_{\mathrm{eff}}(m_{\mathrm{DM}},z)\nonumber \\&\quad =\int ^{m_{\mathrm{DM}}}_0 dE\frac{E}{m_{\mathrm{DM}}}\left[ 2\frac{dN^i_{e^+}(m_{\mathrm{DM}},E)}{dE}f^{e^+e^-}_{\mathrm{eff}}(E,z)\nonumber \right. \\&\qquad \left. +\frac{dN^i_\gamma (m_{\mathrm{DM}},E)}{dE}f^\gamma _{\mathrm{eff}}(E,z)\right] . \end{aligned}$$

(49)

The detailed calculation of $f_{\mathrm{eff}}(z)$ has been developed in Refs. [188, 189, 206,207,208,209,210,211] and the numerical results available at [212] for all of 28 SM final states based on annihilation spectrum provided by PPPC4DMID package [213]. In our model, N annihilate into muons, taus, and neutrinos according to channels $NN\rightarrow Z^\prime Z^\prime \rightarrow 2\ell ^+2\ell ^-+2\nu _\ell 2\nu _\ell$ and $NN\rightarrow Z^\prime H_0(\rightarrow Z'Z') \rightarrow 3\ell ^+3\ell ^-+3\nu _\ell 3\nu _\ell$, ($\ell =\mu ,\tau$). Notice that both $Z^\prime$ and $H_0$ are on-shell mediators in relative annihilation channels. Corresponding $f_{\mathrm{eff}}(z)$ for our benchmarks can be obtained by applying simple kinematics, and expressed in terms of Eq. (49) as follows

Fig. 10

$f_{\mathrm{eff}}(z)$ curves in the range of $z\in 100-1200$ for benchmarks in Tables 3 (left) and 4 (right)

$$\begin{aligned} f^{Z'Z'}_{\mathrm{eff}} (z)= & {} \sum _{i=\mu ,\tau } f^i_{\mathrm{eff}}(M_N/2,z), \nonumber \\ f^{Z'H_0}_{\mathrm{eff}}(z)= & {} \sum _{i=\mu ,\tau }\left[ \frac{f^i_{\mathrm{eff}}(E_{Z'}/2,z)+2(E_{H_0}/E_{Z'})f^i_{\mathrm{eff}}(E_{H_0}/4)}{1+2(E_{H_0}/E_{Z'})}\right] ,\nonumber \\ \end{aligned}$$

(50)

respectively for two annihilation channels. In above equation, $E_{Z',H_0}=M_N[1\pm (M^2_{Z'}-M^2_{H_0})/4M^2_N]$. $f_{\mathrm{eff}}(z)$ curves for benchmarks in Tables 3 and 4 are displayed in Fig. 10 for redshift $z\in 100{-}1200$. Here we only care about $\mu ,~\tau$ final states since neutrinos have negligible energy deposition. In the later case, neutrino detectors such as IceCube can impose a better constraint, but current limits are still weak, thus do not threaten our model.

Based on these curves, We now calculate CMB limits on annihilation cross section. As is shown in previous studies by using the method of principal component analysis, the CMB constraint is dominated by the behavior of $f_{\mathrm{eff}}(z)$ at $z\sim 600$. This then impose upper limit on $p_{\mathrm{ann}}$ function as [169],

$$\begin{aligned} p_{\mathrm{ann}}=&f_{\mathrm{eff}}(z=600)\frac{\mathrm{BR}_{\mu ,\tau }\langle \sigma v\rangle _{\mathrm{CMB}}}{M_N}<3.4\nonumber \\&\times 10^{-28}~\mathrm{cm}^3/\mathrm{s}/\mathrm{GeV}(95\%~\mathrm{C.L.}), \end{aligned}$$

(51)

where $\mathrm{BR}_{\mu ,\tau }\simeq 1/3$ with neglecting phase space difference. Resulting upper bound for annihilation cross section, $\langle \sigma v\rangle _{\mathrm{CMB}}$, are also listed in Tables 3 and 4. Comparing with $\langle \sigma v\rangle _{BF}$, we found that they have slight conflict. However, considering current measurement allowed a larger local density within uncertainty, which will result in a smaller cross section (by a factor of several times) required by positron excess. In this case, our benchmarks are still comparable with CMB constraint marginally.

6 Conclusion

In light of recent $R_{K^{(*)}}$ anomaly and AMS-02 positron excess, we revise the gauged $U(1)_{L_\mu -L_\tau }$ scotogenic model. This model implement gauged $U(1)_{L_\mu -L_\tau }$ symmetry in the one-loop radiative neutrino mass model, where three right-handed neutrinos $N_{\ell }(\ell =e,\mu ,\tau )$, a scalar doublet $\eta$, and a scalar singlet S with charge $+1$ under $U(1)_{L_\mu -L_\tau }$ are introduced. The $U(1)_{L_\mu -L_\tau }$ symmetry is broken spontaneously by the VEV of S, resulting the massive gauge boson $Z'$. As a complementarity for previous consideration of light $Z'$, we instead mainly consider the case of heavy $Z'$, i.e., $M_{Z'}\gtrsim 10~\mathrm{GeV}$.

Provided the existence of certain vector-like quarks charged under $U(1)_{L_\mu -L_\tau }$, an effective $Z' bs$ coupling could be generated, thus the gauge boson $Z'$ will contribute to the process $b\rightarrow s \mu ^+\mu ^-$. In the scenario of heavy $Z'$, the required Wilson coefficient $C_9^\mu \approx -0.95$ to explain $R_{K^{(*)}}$ anomaly can be acquired with Yukawa coupling $Y_Q=0.122$ and $M_{Q}=10~\mathrm{TeV}$. Meanwhile, constraints on $Z'$ mainly come from neutrino trident production and $B_s$ mixing, which actually require $550~\mathrm{GeV}\lesssim M_{Z'}/g'\lesssim 4~\mathrm{TeV}$. For $Z'\lesssim 50~\mathrm{GeV}$, the search for $Z'$ in the $Z\rightarrow 4\mu$ final states set a more tight constraint than neutrino trident production.

As for fermion DM N, the gauge boson $Z'$ and scalar singlet $H_0$ arising from the gauged $U(1)_{L_\mu -L_\tau }$ provide viable annihilation channels. To illustrate the effects of various annihilation processes, we implement a random scan over certain parameter space under constraints from relic density and direct detection. For $M_{Z'}>M_N$, the $NN\rightarrow Z^{\prime *}\rightarrow \ell ^+\ell ^-,\nu _\ell \nu _\ell (\ell =\mu ,\tau )$ are the dominant channel, while $NN\rightarrow Z'Z', Z'H_0,H_0H_0$ channels become dominant for $M_{Z',H_0}<M_N$. Especially, the $NN\rightarrow Z'Z', Z'H_0(\rightarrow Z'Z')$ channels can account for observed positron excess.

Finally, our combined analysis with previous two part shows that the $R_{K^{(*)}}$ anomaly and AMS-02 positron excess can be explained simultaneously under constraints from neutrino trident production, $B_s$ mixing, neutrino mixing, DM relic density, direct detection as well as various indirect detection (AMS-02 antiproton, Fermi-LAT EGRB and CMB measurements).