#  Scalar dark matter and Muon $g-2$ in a $U(1)_{L_{\mu}-L_{\tau}}$ model

XinXin Qi AiGeng Yang Institute of Theoretical Physics, School of Physics,
Dalian University of Technology,
No.2 Linggong Road, Dalian, Liaoning, 116024, People’s Republic of China  Wei
Liu Department of Applied Physics, Nanjing University of Science and
Technology,
Nanjing 210094, People’s Republic of China  Hao Sun<EMAIL_ADDRESS>Institute of Theoretical Physics, School of Physics, Dalian University of
Technology,
No.2 Linggong Road, Dalian, Liaoning, 116024, People’s Republic of China

###### Abstract

We consider a simple scalar dark matter model within the frame of gauged
$L_{\mu}-L_{\tau}$ symmetry. A gauge boson $Z^{\prime}$ as well as two scalar
fields $S$ and $\Phi$ are introduced to the Standard Model (SM). $S$ and
$\Phi$ are SM singlet but both with $U(1)_{L_{\mu}-L_{\tau}}$ charge. The real
component and imaginary component of $S$ can acquire different masses after
spontaneously symmetry breaking, and the lighter one can play the role of dark
matter which is stabilized by the residual $Z_{2}$ symmetry. A viable
parameter space is considered to discuss the possibility of light dark matter
as well as co-annihilation case, and we present current $(g-2)_{\mu}$ anomaly,
Higgs invisible decay, dark matter relic density as well as direct detection
constriants on the parameter space.

##  I Introduction

Recently, the Fermi-lab reported a new measurement from the Muon $g-2$
experiment Abi et al. (2021); Albahri et al. (2021a, b), the result is:

$\displaystyle a_{\mu}^{\rm exp}=116592061(41)\times 10^{-11},$ (1)

compared with the SM prediction

$\displaystyle a_{\mu}^{\rm SM}=116591810(43)\times 10^{-11}.$ (2)

The observed $4.2\sigma$ discrepancy raises a big challenge to the SM and can
give us new hints to the new physics. Models related with $(g-2)_{\mu}$
anomaly can be found in Kawamura et al. (2020); Athron et al. (2021); Baltz
and Gondolo (2001); Everett et al. (2001); Choudhury et al. (2001); Cheung
(2001); Calmet and Neronov (2002); Xiong and Yang (2001); Cox et al. (2021);
Chakrabarty et al. (2018); Chakrabarty (2020). Among these models, the gauged
$U(1)_{L_{\mu}-L_{\tau}}$ has been discussed for a long time for its
simplicity since anomaly cancellation can be accomplished without introducing
new fermions in such model. Discussion about models for a new gauge boson
($Z^{\prime}$) can be found in Foot et al. (1994); Altmannshofer et al.
(2016a); Foldenauer and Jaeckel (2017); Baek et al. (2001); He et al. (1991a,
b); Bélanger et al. (2015); Heeck and Rodejohann (2011), where the
$(g-2)_{\mu}$ anomaly can be naturally explained by the one loop contribution
of $Z^{\prime}$, and the new gauge boson mass is limited to be not too heavy.
Searches for $Z^{\prime}$ at experiment can give the most stringent limit on
the gauge boson mass and coupling constant, and a viable parameter space with
$Z^{\prime}$ boson mass at MeV scale has been well studied, which also
satisfies current experiment constraints Lees et al. (2016); Mishra et al.
(1991); Kaneta and Shimomura (2017); Chun et al. (2019).

The gauged $U(1)_{L_{\mu}-L_{\tau}}$ symmetry gives possible direction for the
extension of SM, and one can introduce new extra particles with
$U(1)_{L_{\mu}-L_{\tau}}$ charge to the model to explain dark matter problem,
which is also a huge challenge to the SM. According to the astronomical
observation, our universe is not just composed of SM particles, but also dark
matter as well as dark energy Komatsu et al. (2011). Dark matter particles are
assumed to be electrically neutral, colorless and stable on cosmological
scales. To explain the observed relic density of dark matter, thermal Freeze-
out Chiu (1966) and Freeze-in Hall et al. (2010) mechanism have been put
forward. Aside from dark matter annihilating into SM particles, other
processes such as semi-annihilation Bélanger et al. (2014) and co-annihilation
Baker et al. (2015), can also contribute to the relic density of dark matter.
New fermions as dark matter candidate in a gauged $L_{\mu}-L_{\tau}$ model can
be found in Foldenauer (2019); Altmannshofer et al. (2016b), while scalar dark
matter models have been discussed in Biswas et al. (2017); Das et al. (2021),
where the current relic density is obtained via Freeze-in mechanism. Among
these models, dark matter particles are almost stabilized by
$U(1)_{L_{\mu}-L_{\tau}}$ symmetry. Besides these models, one can also
consider a type of so-called darkon DM matter model Cheung et al. (2014); He
et al. (2009). The darkon can play the role of dark matter via its lighter
component after spontaneously symmetry breaking in the case of $U(1)$
extension of the SM Chiang et al. (2013); Farzan and Akbarieh (2012, 2014).
Particularly, one can have co-annihilation contribution when the real
component and imaginary component of darkon are nearly degenerate.

In this article, we consider a simple scalar dark matter model within the
frame of the gauged $L_{\mu}-L_{\tau}$ symmetry. We introduce a new gauge
boson $Z^{\prime}$ and two scalar fields $S$ and $\Phi$ to the SM, $S$ and
$\Phi$ are singlet in the SM but carry $U(1)_{L_{\mu}-L_{\tau}}$ charge. The
scalar $\Phi$ will spontaneously breaking so that the new gauge boson
$Z^{\prime}$ acquires mass, while the lighter component of $S$ will play the
role of dark matter in our model after spontaneously breaking, which will be
stabilized with the residual $Z_{2}$ symmetry. We focus on the $(g-2)_{\mu}$
anomaly problem and scalar dark matter in this work, and discussion about
neutrino mass problem in a similar gauged ${L_{\mu}-L_{\tau}}$ model can be
found in Biswas et al. (2016a), where the scalar dark matter in this model is
stabilized by $U(1)_{L_{\mu}-L_{\tau}}$ symmetry. We consider the possiblity
of light dark matter as well as co-annihilation case in our model and we show
the Higgs invisble decay, relic density constraint as well as direct detection
constraint on a viable parameter space.

This article is arranged as followed: We give the description of the scalar
dark matter model in section II. We consider the $(g-2)_{\mu}$ anomaly
constraint on the parameter space in section III. We discuss Higgs invisible
decay, relic density as well as direct detection constraint on the chosen
parameter space separately in IV.1, IV.2 and IV.3. We give a summary in the
last section V.

##  II model description

We discuss a scalar dark matter model based on the $U(1)_{L_{\mu}-L_{\tau}}$
extension of the SM. A new gauge boson $Z^{\prime}$ as well as two SM singlet
scalar $\Phi$ and $S$ are introduced to the standard model. The field $H$
gives the non-zero vaccum expection value (vev) like in the SM, while $\Phi$
develops a vev to break the $U(1)_{L_{\mu}-L_{\tau}}$ symmetry. The field $S$
has zero vev and plays the role of scalar dark matter via its lighter
component. $Z^{\prime}$ will acquire mass after symmetry breaking
spontaneously. We can assume $S$ and $\Phi$ carry opposite
$U(1)_{L_{\mu}-L_{\tau}}$ charge so that we have
$\lambda_{ds}(S^{2}\Phi^{2}+{S^{\dagger}}^{2}{\Phi^{\dagger}}^{2})$ in the
lagrangian, which is essential because it triggers the
$U(1)_{L_{\mu}-L_{\tau}}\to Z_{2}$ spontaneously breakdown as soon as $\Phi$
acquires a vev.

The additional part of the fermion lagrangian is given by:

$\displaystyle{\cal L}_{\rm fermion}$
$\displaystyle=g_{p}Z^{\prime\sigma}(\bar{\mu}\gamma^{\sigma}\mu-\bar{\tau}\gamma^{\sigma}\tau$
$\displaystyle+\bar{v}_{\mu}\gamma^{\sigma}P_{L}v_{\mu}-\bar{v}_{\tau}\gamma^{\sigma}P_{L}v_{\tau}).$
(3)

The scalar part with dark matter is given by:

$\displaystyle{\cal L}_{\rm
scalar}=|D_{\mu}\Phi|^{2}+|D_{\mu}H|^{2}+|D_{\mu}S|^{2}-{\cal V}(H,\Phi,S)$
(4)

with

$\displaystyle D_{\mu}\Phi=\partial_{\mu}\Phi-
ig_{p}q_{1}A^{\prime}_{\mu}\Phi,$ $\displaystyle
D_{\mu}H=\partial_{\mu}H-ig_{L}W_{\mu}H-\frac{i}{2}g_{Y}A_{\mu}H,$ (5)
$\displaystyle D_{\mu}S=\partial_{\mu}S-ig_{p}q_{2}A^{\prime}_{\mu}S,$

where $q_{1}$ and $q_{2}$ are the $U(1)_{L_{\mu}-L_{\tau}}$ charge of $\Phi$
and $S$, $g_{p}$ is the $U(1)_{L_{\mu}-L_{\tau}}$ coupling constant. For
simplicity, we assume $q_{1}=2$ and $q_{2}=-2$ as we considered above. The
potential term ${\cal V}(H,\Phi,S)$ is given by:

$\displaystyle{\cal V}(H,\Phi,S)$
$\displaystyle=\mu|H|^{2}+\lambda_{H}|H|^{4}+\mu_{1}|\Phi|^{2}+\lambda_{p}|\Phi|^{4}$
$\displaystyle-\frac{1}{2}m_{0}^{2}|S|^{2}+\frac{\lambda_{s}}{4}|S|^{4}+\lambda_{Hp}|H|^{2}|\Phi|^{2}$
$\displaystyle+\lambda_{Hs}|H|^{2}|S|^{2}+\lambda_{SP}|S|^{2}|\Phi|^{2}$
$\displaystyle+\frac{\lambda_{ds}}{2}(S^{2}\Phi^{2}+{S^{\dagger}}^{2}{\Phi^{\dagger}}^{2}).$
(6)

where $\mu<0$ and $\mu_{1}<0$ and $m_{0}^{2}>0$. We assume $H$ and $\Phi$
spontaneously breaking, and we have $H$ and $\Phi$ in unitary gauge form with:

$\displaystyle H=\left(\begin{array}[]{c}0\\\
\frac{v+h}{\sqrt{2}}\end{array}\right),~{}~{}~{}\Phi=\frac{v_{b}+\phi}{\sqrt{2}}$
(9)

where $v$ and $v_{b}$ are the vevs and we assume $v$ is the SM one equal 246
GeV. What’s more, $v_{b}=M_{Z_{p}}/2g_{p}$, where $M_{Z_{p}}$ is the new gauge
boson $Z^{\prime}$ mass. Furthermore, we can write $S=S_{R}+iS_{I}$ in the
form of real component and imaginary component, and we have:

$\displaystyle m_{R}^{2}$ $\displaystyle=$ $\displaystyle
m_{0}^{2}+\lambda_{SP}v_{b}^{2}+\lambda_{Hs}v^{2}+\lambda_{ds}v_{b}^{2}$ (10)
$\displaystyle m_{I}^{2}$ $\displaystyle=$ $\displaystyle
m_{0}^{2}+\lambda_{SP}v_{b}^{2}+\lambda_{Hs}v^{2}-\lambda_{ds}v_{b}^{2}$ (11)

where $m_{R}^{2}$, $m_{I}^{2}$ correspond to the squared mass of $S_{R}$ and
$S_{I}$. The mass difference between $S_{R}$ and $S_{I}$ is determined by the
sign of $\lambda_{ds}$. We can take $\lambda_{ds}>0$ so that the lighter
particle $S_{I}$ is the weakly interacting massive particle (WIMP) dark
matter. The result will be the same if we take $\lambda_{ds}<0$ while $S_{R}$
acts as the dark matter.

The mass matrix related with two Higgs is given by the following, after
spontaneously breaking:

$\displaystyle{\cal
M}=\left(\begin{array}[]{cc}2\lambda_{H}v^{2}&\lambda_{Hp}vv_{b}\\\
\lambda_{Hp}vv_{b}&2\lambda_{p}v^{2}_{b}\\\ \end{array}\right).$ (14)

The mass matrix eigenvalue can be analytically expressed, and the result is
given by:

$\displaystyle m_{1}^{2}$ $\displaystyle=$
$\displaystyle\lambda_{H}v^{2}+\lambda_{p}v_{b}^{2}-\sqrt{(\lambda_{H}v^{2}-\lambda_{p}v_{b}^{2})^{2}+\lambda_{Hp}^{2}v_{b}^{2}v^{2}}~{}~{}~{}~{}$
(15) $\displaystyle m_{2}^{2}$ $\displaystyle=$
$\displaystyle\lambda_{H}v^{2}+\lambda_{p}v_{b}^{2}+\sqrt{(\lambda_{H}v^{2}-\lambda_{p}v_{b}^{2})^{2}+\lambda_{Hp}^{2}v_{b}^{2}v^{2}}.~{}~{}~{}~{}~{}~{}$
(16)

The gauge eigenstate and mass eignestate is related with a mixing angle
$\theta$ which can be determined

$\displaystyle\sin
2\theta=\frac{\lambda_{Hp}v_{b}v}{\sqrt{(\lambda_{H}v^{2}-\lambda_{p}v^{2}_{b})^{2}+\lambda_{Hp}^{2}v_{b}^{2}v^{2}}}.$
(17)

We consider the lighter Higgs is just the SM Higgs observed at the LHC with
$m_{1}=125$GeV and $m_{2}$ being another Higgs mass in this work. We use
$h_{1}$ to represent the SM Higgs and $h_{2}$ to represent another Higgs for
convenience. The relation between the mass eigenstates and mass eigenstates
can be given by:

$\left(\begin{array}[]{c}h_{1}\\\ h_{2}\\\
\end{array}\right)=\left(\begin{array}[]{cc}\cos\theta&-\sin\theta\\\
\sin\theta&\cos\theta\\\ \end{array}\right)\left(\begin{array}[]{c}h\\\
\phi\\\ \end{array}\right)$ (18)

and we can express the couplings with two Higgs mass and mixing angle as
followed:

$\displaystyle\lambda_{H}=\frac{m_{1}^{2}}{4v^{2}}(1+\cos
2\theta)+\frac{m_{2}^{2}}{4v^{2}}(1-\cos 2\theta)$ (19)
$\displaystyle\lambda_{p}=\frac{m_{1}^{2}}{4v_{b}^{2}}(1-\cos
2\theta)+\frac{m_{2}^{2}}{4v_{b}^{2}}(1+\cos 2\theta)$ (20)
$\displaystyle\lambda_{Hp}=\sin 2\theta\frac{m_{2}^{2}-m_{1}^{2}}{2vv_{b}}.$
(21)

To make sure the model perturbative, contribution from loop correction should
be smaller than tree level values, such constraint can be ensured when
Chakrabarty et al. (2015)

$\displaystyle|\lambda_{H}|<4\pi,\ |\lambda_{Hp}|<4\pi,\
|\lambda_{p}|<4\pi,|\lambda_{Hs}|<4\pi,$ $\displaystyle|\lambda_{ds}|<4\pi,\
|\lambda_{SP}|<4\pi.$ (22)

What’s more, to obtain a stable vaccum, the quartic couplings apprearing in
the lagrangian should be constrained, we consider the necessary and sufficient
conditions as followed Kannike (2016); Biswas et al. (2016b):

$\displaystyle\lambda_{H}\geq 0,\lambda_{s}\geq 0,\lambda_{p}\geq
0,\lambda_{Hs}\geq-2\sqrt{\lambda_{H}\lambda_{s}},$
$\displaystyle\lambda_{Hp}\geq-2\sqrt{\lambda_{H}\lambda_{p}},\
\lambda_{SP}\geq-2\sqrt{\lambda_{P}\lambda_{s}},$
$\displaystyle\sqrt{\lambda_{Hs}+2\sqrt{\lambda_{H}\lambda_{s}}}\sqrt{\lambda_{Hp}+2\sqrt{\lambda_{H}\lambda_{p}}}\sqrt{\lambda_{SP}+2\sqrt{\lambda_{p}\lambda_{s}}}$
$\displaystyle+2\sqrt{\lambda_{H}\lambda_{s}\lambda_{p}}+\lambda_{Hp}\sqrt{\lambda_{s}}+\lambda_{Hs}\sqrt{\lambda_{p}}$
$\displaystyle+\lambda_{SP}\sqrt{\lambda_{H}}\geq 0.$ (23)

In this work, we choose the following parameters as the inputs with

$\displaystyle m_{2},\ M_{Z_{p}},\ g_{p},\ \lambda_{s},\ \lambda_{ds},\
\lambda_{SP},\ \lambda_{Hs},\ \sin\theta,\ m_{0}.$ (24)

##  III Parameter constraints on $Z^{\prime}$

Figure 1: One-loop contribution of $Z^{\prime}$ to $(g-2)_{\mu}$ anomaly.

Before discussion, we should consider the $(g-2)_{\mu}$ constraints on
$M_{Z_{p}}-g_{p}$ plane, since one prospect to introduce
$U(1)_{L_{\mu}-L_{\tau}}$ symmetry is to explain the $(g-2)_{\mu}$ anomaly.
Contribution of the new gauge boson $Z^{\prime}$ to the muon anomalous
magnetic momentum is shown in Figure 1. The analytical expression for $\Delta
a_{\mu}$ is

$\displaystyle\Delta
a_{\mu}=\frac{g_{p}^{2}m_{\mu}^{2}}{8\pi^{2}}\int_{0}^{1}\frac{2x^{2}(1-x)}{x^{2}m_{\mu}^{2}+(1-x)M_{Z_{p}}^{2}}dx$
(25)

where $m_{\mu}$ is the muon mass.

Searches for new gauge bosons at experiment has already given the most
stringent constraint on the parameter space of $M_{Z_{p}}-g_{p}$. At
colliders, the $Z^{\prime}$ gauge boson can be produced via
$e^{+}e^{-}\to\mu^{+}\mu^{-}Z^{\prime}$ with the subsequent decay of
$Z^{\prime}\to\mu\mu$. Such search can be found in Babar Lees et al. (2016)
and give us the possible bound on $M_{Z_{p}}-g_{p}$ plane. What’s more,
neutrino expriments give us new clues to constrain the parameter space. The
Borexino data related with the scattering of low energy solar neutrinos
Agostini et al. (2019); Abdullah et al. (2018) can provide the most stringent
constraint on the low $M_{Z_{p}}$ and low $g_{p}$ region. Another constraint
is from CCFR collaboration Altmannshofer et al. (2014), obtained via the
neutrino trident production which is related with the process $\nu N\to\nu
N\mu^{+}\mu^{-}$, where a muon neutrino scattered off of a nucleus producing a
$\mu^{+}\mu^{-}$ pair. Such process will be enhanced due to the existence of
$Z^{\prime}$ compared with SM case, which gives strong bounds on the possible
$Z^{\prime}$ contribution and constrain the $M_{Z_{p}}-g_{p}$ parameter space
accordingly.

Figure 2:  The allowed region that can explain the $(g-2)_{\mu}$ anomaly and
the exclued region on $M_{Z_{p}}-g_{p}$ from different experiments: the green
(orange, blue, purple) region is excluded by the Borexino (BABAR, LHC, CCFC).
The red region is the parameter satisfying $(g-2)_{\mu}$ anomaly.

In our paper, we focus on the low $M_{Z_{p}}$ region. The combined experiment
as well as $(g-2)_{\mu}$ constraints on the $M_{Z_{p}}-g_{p}$ is given in
Figure 2. According to Figure 2, CCFR gives the most stringent constraint on
the parameter space we interest at low $M_{Z_{p}}$ region. The Borexino limits
become relavant at smaller $M_{Z_{p}}$ and $g_{p}$. For
$M_{Z_{p}}\subseteq[0.2,5]$ GeV, BABAR data plays an important role to
constrain $M_{Z_{p}}$ and $g_{p}$.

Figure 3:  $M_{Z_{p}}-g_{p}$ plane to explain $(g-2)_{\mu}$ anomaly with
$M_{Z_{p}}\subseteq[0.01,0.25]$ GeV. The red region is Log result of Eq.(25)
and the two yellow regions are the latest experiment result of $(g-2)_{\mu}$.
The overlap part of the three regions is the viable parameter space to explain
$(g-2)_{\mu}$ anomaly within experiment constraints.

To sum up, most of the parameter space to explain $(g-2)_{\mu}$ has been
excluded by these experiments but $M_{Z_{p}}\subseteq[0.01,0.25]$ GeV and
$g_{p}$ at $10^{-3}$ level. The viable parameter region that can explain the
$(g-2)_{\mu}$ anomaly in the $M_{Z_{p}}-g_{p}$ plane is shown in Figure 3
accoding to the latest experiment value $\Delta a_{\mu}=251(59)\times
10^{-11}$ of where the red plane is the result of $\rm Log\Delta a_{\mu}$ and
the two yellow palnes are the low bound and top bound on $\rm Log\Delta
a_{\mu}$ from the experiment. Intersection region of the three planes is the
viable parameter space of $M_{Z_{p}}-g_{p}$ explaining $(g-2)_{\mu}$ anomaly,
and we can give some pairs of $(M_{Z_{p}},g_{p})$ with $(0.01\rm
GeV,4.18\times 10^{-4})$, $(0.25\rm GeV,1.117\times 10^{-3})$.

##  IV Phenomenological study

###  IV.1 Higgs invisible decay

In our model, we have a new gauge boson $Z^{\prime}$ and two scalar fields
$S_{I}$ and $S_{R}$, We assume all these particles masses are smaller than the
SM Higgs mass with $M_{Z_{p}}<1/2m_{1}$, $m_{I}<1/2m_{1}$ and $m_{R}<1/2m_{1}$
so that SM Higgs can decay to these particles which will contribute to the
invisible Higgs decay width. Current constraint on such invisible branching
fraction is ${\rm Br}(h\to inv)<0.24$ accroding to the observations at the LHC
for the SM Higgs Tanabashi et al. (2018), which means

$\displaystyle\frac{\Gamma(h\to inv)}{\Gamma(h\to inv)+\Gamma(h\to{\rm
SM})}<0.24.$ (26)

The decay widths related with Higgs invisible decay channel in our model are
given as followed:

$\displaystyle\Gamma_{h_{1}\to S_{I}S_{I}}$ $\displaystyle=$
$\displaystyle\sqrt{m_{1}^{2}(m_{1}^{2}-4m_{I}^{2})}/(8\pi m_{1}^{3})$
$\displaystyle\times$
$\displaystyle(\cos\theta\lambda_{Hs}v+\lambda_{ds}\sin\theta
v_{b}-\lambda_{SP}\sin\theta v_{b})^{2}$ $\displaystyle\Gamma_{h_{1}\to
S_{R}S_{R}}$ $\displaystyle=$
$\displaystyle\sqrt{m_{1}^{2}(m_{1}^{2}-4m_{R}^{2})}/(8\pi m_{1}^{3})$
$\displaystyle\times$
$\displaystyle(\cos\theta\lambda_{Hs}v-\lambda_{ds}\sin\theta
v_{b}-\lambda_{SP}\sin\theta v_{b})^{2}$ $\displaystyle\Gamma_{h_{1}\to
Z^{\prime}Z^{\prime}}$ $\displaystyle=$ $\displaystyle
g_{p}^{4}\sqrt{m_{1}^{2}(m_{1}^{2}-4M_{Z_{p}}^{2})}/(2M_{Z_{p}}^{4}\pi
m_{1}^{3})$ (27) $\displaystyle\times$
$\displaystyle(m_{1}^{4}-4m_{1}^{2}M_{Z_{p}}^{2}+12M_{Z_{p}}^{4})\sin^{2}\theta
v_{b}^{2}.$

Figure 4: Contourplot of Higgs invisible decay width where the X-axis is
$m_{0}$ and Y-axis is $\lambda_{ds}$ when $\lambda_{Hs}=0.005$,
$\lambda_{SP}=0.001$.

In this part, we consider the Higgs invisible decay constraint on the chosen
parameter space. As we can see in the following discussion, $\lambda_{ds}$ is
limited stringetly under relic density constraint, while $\lambda_{Hs}$ and
$\lambda_{SP}$ are more flexible. According to Eq.(26), the contribution of
$\lambda_{Hs}$ to Higgs invisble decay width is positive while $\lambda_{SP}$
is negative. For simplicity, we can set $\lambda_{Hs}$ and $\lambda_{SP}$ to
be some special value and give the result of the Higgs invisible decay width
of our model to estimate the chosen parameter space. In Figure 4, we fixed
$\lambda_{Hs}=0.005$ and $\lambda_{SP}=0.001$ and set
$\lambda_{ds}\subseteq[1\times 10^{-5},0.01]$ and $m_{0}\subseteq[0,60]GeV$.
According to Figure 4, the Higgs invisible decay width is much lower than
current constraint in the case of $\lambda_{Hs}=0.005$ and
$\lambda_{SP}=0.001$. As we discussed above, we can set $\lambda_{Hs}<0.005$
and $\lambda_{SP}>0.001$ to get a viable parameter space satisfying Higgs
invisible decay width constraint.

###  IV.2 Relic density

In this part, we discuss the dark matter phenomenology in the model. The
expression of relic abundance of dark matter can be given as followed:

$\displaystyle\Omega h^{2}=\frac{1.07\times 10^{9}{\rm
GeV}^{-1}}{{g*}^{1/2}M_{Pl}}\frac{1}{J(x_{f})}$ (28)

where $g*$ is the total number of effective relativistic degrees of freedom,
$M_{Pl}=1.22\times 10^{19}$ GeV is the Planck mass, and $J(x_{f})$ is given
by:

$\displaystyle J(x_{f})=\int_{x_{f}}^{\infty}\frac{dx}{x^{2}}\langle\sigma{\rm
v}\rangle(x).$ (29)

The freeze-out parameter $x_{f}$ in the integral is Kolb and Turner (1981):

$\displaystyle x_{f}=\ln\frac{0.038gM_{Pl}m_{\rm DM}\langle\sigma{\rm
v}\rangle(x_{f})}{({g*}x_{f})^{1/2}}$ (30)

where $m_{\rm DM}$ is the dark matter mass. According to current experiment,
the dark matter relic density is Aghanim et al. (2020):

$\displaystyle\Omega h^{2}=0.1198\pm 0.0012.$ (31)

Scalar dark matter in our model is stabilized by $Z_{2}$ symmetry after
$U(1)_{L_{\mu}-L_{\tau}}$ symmetry spontaneously breaking. It is worth
stressed that another scalar $S_{R}$ can also play role of dark matter when
$S_{I}$ and $S_{R}$ are nearly degenerate. The dark matter relic density will
not just be determined by $S_{I}$ but also $S_{R}$, which is so-called co-
annihilation case Griest and Seckel (1991). In the case we considered, the
number density of $S_{R}$ will track $S_{I}$ number density during freeze-out,
when the relative mass splitting $\Delta$ is small compared to the freeze-out
temperature, which is defined by:

$\displaystyle\Delta\equiv\frac{m_{R}-m_{I}}{m_{I}}.$ (32)

Concretely speaking, the relic density is calculated by solving the Boltzman
equation, which depends on the dimensionless variable $x=m_{\rm DM}/T$.
Freeze-out occurs at $x_{F}=m_{\rm DM}/T_{F}\approx 20\sim 30$ for cold, non-
relativistic dark matter, where $T_{F}$ is the freeze-out temperature. For
$\delta m=m_{R}-m_{I}$ much larger than $T_{F}$, we have just dark matter
$S_{I}$ annihilations freeze-out. However, $S_{R}$ will be thermally
accessible when $\delta m\approx T_{F}$. Hence, the relative mass splitting
$\Delta$ can give the upper bound that $\Delta\sim x_{F}^{-1}\sim 0.03-0.05$
when we consider the co-annihilation case. A more systematic estimate for the
contribution of $\Delta$ to relic density can be found in Baker et al. (2015).
$\Delta$ in our model can be given as followed:

$\displaystyle\Delta=\frac{m^{2}_{R}-m^{2}_{I}}{m_{I}(m_{I}+m_{R})}=\frac{2\lambda_{ds}v^{2}_{b}}{m^{2}_{I}+m_{I}m_{R}}.$
(33)

According to Eq.(33), co-annihilation becomes more significant in the model as
long as $\lambda_{ds}\ll m_{I}$, where $S_{I}$ and $S_{R}$ are kept in
equilibrium via the interactions $S_{I}S_{I}\leftrightarrow S_{R}S_{R}$, while
in the light dark matter case the main contribution of relic density arising
from the annihilation of $S_{I}$. Before we consider the relic density
numberically, we should stress that co-annihilation process can be dominated
during the evolution of relic denstiy in the heavy dark matter case, but we
want to consider a viable region where annihilation and co-annihilation can be
both involved which means the dark matter mass we considered should not be too
heavy. We choose a viable parameter space with $0\leqslant m_{0}\leqslant 60$
GeV and dark matter particle mass will be constrained to be a few GeV to about
70 GeV. Annihilation process and co-annihilation will be both involved within
the chosen parameter space as we can see in the following discussion.

The Feynman diagrams relavant for the dark matter production are given in
Figure 5.

Figure 5:  Feynman diagrams related with dark matter relic density.

According to Figure 5, scalar dark matter can annihilate into SM particles
with Higgs-mediated interactions. In the case of $2m_{\rm DM}>M_{Z_{p}}$ as we
have assumed, dark matter can annihilate into a pair of $Z^{\prime}$ via
t-channel as shown in Figure 5(e). Vertices related with these Feynman
diagrams are given in Table 1 for $S_{I}$ annihilation. As we have discussed
above, $S_{R}$ can also plays the role of dark matter in the case of co-
annihilation, the relavant vertices are also given in Table 1.

Coupling | Vertex Factor
---|---
${\cal C}_{S_{I(R)}S_{I(R)}h_{1}}$ | $-i[2v\lambda_{Hs}\cos\theta\pm 4v_{b}\lambda_{ds}\sin\theta\mp 2v_{b}\lambda_{SP}\sin\theta]$
${\cal C}_{S_{I(R)}S_{I(R)}h_{2}}$ | $-i[2v\lambda_{Hs}\sin\theta\mp 4v_{b}\lambda_{ds}\cos\theta\pm 2v_{b}\lambda_{SP}\cos\theta]$
${\cal C}_{S_{I(R)}S_{I(R)}h_{1}h_{1}}$ | $-i[2\lambda_{Hs}\cos^{2}\theta\mp 4\lambda_{ds}\sin^{2}\theta+2\lambda_{SP}\sin^{2}\theta]$
${\cal C}_{S_{I(R)}S_{I(R)}h_{2}h_{2}}$ | $-i[2\lambda_{SP}\cos^{2}\theta\mp 4\lambda_{ds}\cos^{2}\theta+2\lambda_{Hs}\sin^{2}\theta]$
${\cal C}_{S_{I(R)}S_{I(R)}h_{1}h_{2}}$ | $-i\cos\theta\sin\theta[2\lambda_{Hs}\pm 4\lambda_{ds}-2\lambda_{SP}]$
${\cal C}_{S_{I}S_{R}Z^{\prime}}$ | $\frac{2}{3}g_{p}(p_{1}-p_{2})^{\mu}$
${\cal C}_{S_{I(R)}S_{I(R)}Z^{\prime}Z^{\prime}}$ | $i16g_{p}^{2}$
Table 1:  Vertices related with $S_{I(R)}$ annihilation.

In this work, we use Feynrules Alloul et al. (2014) to generate implemented
code, micrOmegas Bélanger et al. (2018) to calculate the relic density and
t3ps Maurer (2016) to scan the parameter space. For simplicity, we assume
$\lambda_{s}=0$ since the self interaction part of dark matter makes no
contribution to the relic density. Moreover, searches for the exotic Higgs at
the LHC give upper limits of the mixing angle with $|\sin\theta|\leqslant
0.2\sim 0.6$ depending on the heavy Higgs mass Accomando et al. (2016); Robens
and Stefaniak (2015); López-Val and Robens (2014). In addition, as we
discussed above, $M_{Z_{p}}$ and $g_{p}$ are limited with current experiments.
To sum up, there are 8 independent parameters in the model, and we take three
Scenarios to estimate these parameters seprately for simplicity. In $Scenario\
A$, we consider $(g-2)_{\mu}$ anomaly as well as dark matter relic density
constraint on the $M_{Z_{p}}-g_{p}$. In $Scenario\ B$, we discuss contribution
of other parameters on the dark matter relic density. In $Scenario\ C$, we
focus on the relic density constraint on the couplings $\lambda_{Hs}$,
$\lambda_{SP}$ and $\lambda_{ds}$.

#### Scenario A

In this part, we focus on $(g-2)_{\mu}$ anomaly as well as dark matter relic
density constraint on the $M_{Z_{p}}-g_{p}$, the inputs are set by Table 2.

Parameter | value for inputs
---|---
$m_{2}$ |  $[0.2,2]$ TeV
$\sin\theta$ | 0.001
$\lambda_{Hs}$ | $[10^{-5},0.1]$
$\lambda_{SP}$ | $[10^{-5},0.1]$
$\lambda_{ds}$ | $[10^{-5},0.1]$
$M_{Z_{p}}$ |  $[0.01,0.25]$ GeV
$g_{p}$ | $[4\times 10^{-4},0.002]$
$m_{0}$ |  $[0,60]$ GeV
Table 2:  Values for the parameters as input to calculate dark matter relic
density.

According to Figure 6, we give the allowed region to satisfy relic density
constriant (blue dots) and $(g-2)_{\mu}$ anomaly (red region) with
$M_{Z_{p}}\subseteq[0.01,0.25]$ GeV and $g_{p}\subseteq[4\times
10^{-4},0.002]$ where we have set $\sin\theta=0.01$, $m_{2}\subseteq[0.2,2]$
TeV and $\lambda_{Hs,SP,ds}\subseteq[10^{-5},0.1]$. Since we have chosen
$m_{2}$ and $v_{b}=M_{Z_{p}}/2g_{p}$ as inputs, these parameters should be
limited by perturbative constraint as well as vacuum stability constriant,
which means $v_{b}=M_{Z_{p}}/2g_{p}$ should not be too small. As we can see
from Figure 6, most of the blue points fall in the lower-right region, and the
upper-left region is excluded within the chosen parameter space.

Figure 6:  Allowed $M_{Z_{p}}-g_{p}$ region to satisfy relic density
constriant [blue dots] and $(g-2)_{\mu}$ anomaly [red region] with
$M_{Z_{p}}\subseteq[0.01,0.25]$ GeV and $g_{p}\subseteq[4\times
10^{-4},0.002]$ where we have set $\sin\theta=0.01$, $m_{2}\subseteq[0.2,2]$
TeV and $\lambda_{Hs,SP,ds}\subseteq[10^{-5},0.1]$, $m_{0}\subseteq[0,60]$
GeV.

#### Scenario B

Figure 7:  The value of Relic density as function of $m_{0}$. The purple solid
curve is the observed relic density at experiment Aghanim et al. (2020).
Curves with different color corresponding to different parameter varies when
$M_{Z_{p}}=0.2$ GeV, $g_{p}=0.001$ [first figure] and $M_{Z_{p}}=0.1$ GeV,
$g_{p}=0.0007$ [second figure].

In this part, we discuss contribution of other parameters to the dark matter
relic density. The results of relic density are shown in Figure 7, where
curves with different color corresponding with one of the parameter varies. We
set $0\leq m_{0}\leq 60\rm GeV$ and benchmark value for the parameters are
given in Table 3. We fix $M_{Z_{p}}=0.2$ GeV and $g_{p}=0.001$ in Figure
7[first], while in Figure 7[second] we fix $M_{Z_{p}}=0.1$ GeV and
$g_{p}=0.0007$.

Parameter | benchmark value for inputs
---|---
$m_{2}$ | 300 GeV
$\sin\theta$ | 0.001
$\lambda_{Hs}$ | $1\times 10^{-4}$
$\lambda_{SP}$ | $1\times 10^{-4}$
$\lambda_{ds}$ | $1\times 10^{-4}$
$M_{Z_{p}}$ | 0.2 GeV
$g_{p}$ | 0.001
Table 3: Benchmark value for the parameters as input to calculate dark matter
relic density.

According to Figure 7[first], the purple solid curve is the observed relic
density at experiments. Curves corresponding to $m_{2}=600$ GeV and
$\sin\theta=0.03$ almost coincide with the benchmark curve, this is due to
dark matter annihilation channels related with Higgs make little contribution
to relic denisty in these cases. Correspondingly, for the red line with
$\lambda_{Hs}=0.01$, we have a resonance region at about $2m_{\rm DM}=m_{1}$,
where the relic density drops sharply, and intersect with the relic density
constraint curve. For $\lambda_{ds}=0.001$, see the blue dashed curve, the
dark matter mass can be negative when $m_{0}$ takes small value, so that the
start point of the curve is not $m_{0}=0$. For Figure 7[second] with
$M_{Z_{p}}=0.1$ GeV and $g_{p}=0.0007$, the curves are almost the same with
those in Figure 7[first], and one of the typical difference is that the blue
dashed curve corresponding to $\lambda_{ds}=0.001$ starts at $m_{0}=0$ and
intersect with the relic density constraint curve. According to these figures,
we can reduce the inputs to $\lambda_{ds}$, $\lambda_{SP}$, $\lambda_{Hs}$,
$m_{0}$, $M_{Z_{p}}$ and $g_{p}$, since contribution of $m_{2}$ and
$\sin\theta$ to relic density can be small.

#### Scenario C

Parameter | value for inputs
---|---
$m_{2}$ | 300 GeV
$\sin\theta$ | 0.001
$M_{Z_{p}}$ | 0.2 GeV
$g_{p}$ | 0.001
$m_{0}$ |  $[0,60]$ GeV
$\lambda_{Hs}$ | $[10^{-5},0.01]$
$\lambda_{SP}$ | $[10^{-5},0.01]$
$\lambda_{ds}$ | $[10^{-5},0.01]$
Table 4:  values for the parameters as input to scan.

In this part, we scan the parameter space to study the relic density
constraint on the couplings. For simplicity, we fix $M_{Z_{p}}$, $g_{p}$,
$\sin\theta$, $m_{2}$ and focus on $\lambda_{SP}$, $\lambda_{ds}$ as well as
$\lambda_{Hs}$. We set these parameters as in Table 4. To avoid $m_{\rm DM}$
taking too large value, we have set the couplings $\lambda_{Hs,SP,ds}$ to be
smaller than 0.01. The results are given in Figure 8.

Figure 8:  Scatter points satisfying relic density constraint with the X-axis
being $S_{I}$ mass. Y-axis being $S_{R}$ mass in the first figure. The scan
results of $\lambda_{ds}-m_{I}$, $\lambda_{Hs}-m_{I}$ and $\lambda_{SP}-m_{I}$
are in the second, third and fourth figure respectively.

According to the first picture of Figure 8, the possible dark matter $S_{I}$
mass ranges from a few GeV to about $m_{1}/2$ as we discussed above, for
$m_{I}$ bigger than about 10 GeV, we have $m_{R}\approx m_{I}$ and
$\Delta\approx 0$ which means co-annihilation process can be domainate in the
evolution of relic density where both $S_{R}$ and $S_{I}$ play the role of
dark matter. For $S_{I}$ takes smaller value, $\Delta$ can be much large and
only $S_{I}$ plays the role of dark matter, but such region is contrained
stringently. In addition, we give the scan result of $\lambda_{ds}-m_{I}$,
$\lambda_{Hs}-m_{I}$ and $\lambda_{SP}-m_{I}$ in the other figures in Figure 8
seperately. As we can see from the second picture, the allowed value of
$\lambda_{ds}$ is constrained within a narrow region when $m_{I}$ smaller than
about 50 GeV, increasing with $m_{I}$ increases, such result seemingly
contraries to Eq.(11), where $m_{I}$ decreases with $\lambda_{ds}$ increase,
this means $\lambda_{Hs}$ and $\lambda_{SP}$ play more important role
determining dark matter mass under the relic density constraint. For
$m_{I}>50$ GeV, $\lambda_{ds}$ can take value within the whole range of
$[1\times 10^{-5},0.0084]$. For $\lambda_{Hs}$ and $\lambda_{SP}$ in the last
two picture, both can take value from $1\times 10^{-5}$ to $0.01$, and a
certain number of the points fall in the region $[0,001,0.01]$, with a few
points obtained at low mass region.

###  IV.3 Direct detection

Experiments related with dark matter direct detection such as CDMS II Ahmed et
al. (2010), XENON100 Aprile et al. (2017), XENON1T Aprile et al. (2016), LUX
Akerib et al. (2014) have been searching for the signal of the interaction of
DM with nucleon. In our model, quarks do not couple with $Z^{\prime}$ but only
with two Higgs particles and these interactions can be concluded by the
scattering of dark matter particle off a SM fermion $f$ via the t-channel
exchange of the two Higgs particles, which is similar with the two singlet
scalar model Abada et al. (2011); Basak et al. (2021). The effective
lagrangian for dark matter-quark elastic scattering can be given by :

$\displaystyle{\cal L}_{q,eff}=-\frac{m_{q}}{2v}(\frac{{\cal
C}_{h_{1}SS}\cos\theta}{m_{1}^{2}}+\frac{{\cal
C}_{h_{2}SS}\sin\theta}{m_{2}^{2}})SS\bar{q}q,$ (34)

and effective lagrangian related with $S_{I}$ and $S_{R}$ can be given by:

$\displaystyle{\cal
L}_{q,eff}=\sum_{S=S_{R},S_{I}}-\frac{m_{q}}{2v}(\frac{{\cal
C}_{h_{1}SS}\cos\theta}{m_{1}^{2}}+\frac{{\cal
C}_{h_{2}SS}\sin\theta}{m_{2}^{2}})SS\bar{q}q.~{}$ (35)

Furthermore, the effective lagrangian related with DM-nucleon elastic
scattering can be given by:

$\displaystyle{\cal L}_{N,eff}$ $\displaystyle=$
$\displaystyle\sum_{S=S_{R},S_{I}}\frac{m_{N}-\frac{7}{9}m_{B}}{v}$ (36)
$\displaystyle(\frac{{\cal C}_{h_{1}SS}\cos\theta}{m_{1}^{2}}+\frac{{\cal
C}_{h_{2}SS}\sin\theta}{m_{2}^{2}})SS\bar{N}N$

where $m_{N}$ represents the nucleon mass and $m_{B}$ represents the baryon
mass in the chiral limit He et al. (2009). The total cross section for S-N
elastic scattering can be given by:

$\displaystyle\sigma_{SN\rightarrow SN}$ $\displaystyle=$
$\displaystyle\frac{m^{4}_{N}f_{N}^{2}}{4\pi(m_{N}+m_{\rm DM})^{2}}$ (37)
$\displaystyle\times$ $\displaystyle(\frac{{\cal
C}_{h_{1}SS}\cos\theta}{m_{1}^{2}}+\frac{{\cal
C}_{h_{2}SS}\sin\theta}{m_{2}^{2}})^{2},$

where $f_{N}$ is the Higgs-nucleon Formfactor with $f_{N}=0.308(18)$ according
to the phenomenological and lattice-QCD calculations Hoferichter et al. (2017)
and $m_{\rm DM}$ corresponds to both $m_{I}$ and $m_{R}$.

Figure 9:  Spin indepedent cross section as a function of $m_{I}$, the red
points represent result of the chosen parameter space satisfying relic density
constraint, the blue dashed curve represents reslut of XENON1T Aprile et al.
(2016).

We conisder the direct detection constraint on the chosen parameter space of
the model and the result is given in Figure 9, where the plot is drawn as a
function of $m_{I}$. The blue dashed curve in Figure 9 represents reslut of
XENON1T Aprile et al. (2016) and the red points represent result of the chosen
parameter space satisfying relic density constraint. According to Figure 9,
although a certain number of the red points fall above the region of direct
detection constraint, which means these points are excluded by direct
detection constraint. There remain points that mass ranging from a few GeV to
about 60 GeV to both satisfy relic density constraint as well as direct
detection constraint. In other word, while direct detection can give stringent
limit on the parameter space, dark matter mass can be from a few GeV to about
60 GeV within a viable parameter space satisfying relic density constraint in
the model.

##  V Summary

SM has achieved great success for its high accurancy to describle electroweak
and strong interaction. However, there remians problems such as neutrino mass
and dark matter that SM can not explain. In addition, recent results about
muon anomalous magnetic moment $(g-2)_{\mu}$ brings new challenges to the SM.
The $4.2\sigma$ discrepancy between experiment and SM prediction seems to
indicate new physics behind $(g-2)_{\mu}$ anomaly and gives possible hints to
the Beyond SM physics. Among these models, a gauged $U(1)_{L_{\mu}-L{\tau}}$
model stands out for its simplicity since the gauge anomaly cancellation can
be accomplished without introducing extra fermions. Related discussion about
such model and experiment constrain the new guage boson $Z^{\prime}$ mass at
MeV scale, which also satisfying the $(g-2)_{\mu}$ constraint. In addition,
one can also introduce new particles to the ${L_{\mu}-L_{\tau}}$ model so that
dark matter problem as well as neutrino mass problem can both be explained.

In this article, we consider a scalar dark matter model within the frame of a
gauged ${L_{\mu}-L_{\tau}}$ model. We introudce a new gauge boson $Z^{\prime}$
as well as two scalar fields $S$ and $\Phi$ to the SM. The lighter component
of $S$ can play the role of dark matter which is stablized by the residual
$Z_{2}$ symmetry after spontaneously symmetry breaking. In this work, we focus
on dark matter and $(g-2)_{\mu}$ anomaly and ignore the neutrino mass problem.
In the case of $m_{I}<m_{1}/2$ , $m_{R}<m_{1}/2$ and $M_{Z_{p}}<m_{1}/2$, we
have new Higgs invisble channels which should be limited by current experiment
result at LHC. A viable parameter space is considered to discuss the
possibility of light dark matter as well as co-annihilation case. We consider
relic density constraint and find in the case of $m_{I}$ taking a few GeV, we
can have the imaginary part $S_{I}$ of $S$ as light dark matter with $m_{R}$
much larger than $m_{I}$, but the allowed parameter space is constrained
stringently. However, with $m_{I}$ increasing, the allowed value for $m_{I}$
and $m_{R}$ is approximately equal so that we have $\Delta\approx 0$, and we
come to so-called co-annihilation process, where both $S_{R}$ and $S_{I}$ play
the role of dark matter. At low mass region, relic density constraint limits
$\lambda_{ds}$ stringently, which plays a significant role determing $m_{I}$
and $m_{R}$. In addition, direct detection has also been taken into
consideration to constrain the chosen parameter space, and the spin-
independent dark matter nucleon elastic scattering cross section can give
stringent constraint the viable parameter space. As we discussed above, we
have shown this model can explain dark matter and $(g-2)_{\mu}$ anomaly at the
same time in certain parameter space.

###### Acknowledgements.

Hao Sun is supported by the National Natural Science Foundation of China
(Grant No.12075043).

## References

  * Abi et al. (2021) B. Abi et al. (Muon g-2), Phys. Rev. Lett. 126, 141801 (2021), eprint 2104.03281.
  * Albahri et al. (2021a) T. Albahri et al. (Muon g-2), Phys. Rev. A 103, 042208 (2021a), eprint 2104.03201.
  * Albahri et al. (2021b) T. Albahri et al. (Muon g-2), Phys. Rev. D 103, 072002 (2021b), eprint 2104.03247.
  * Kawamura et al. (2020) J. Kawamura, S. Okawa, and Y. Omura, JHEP 08, 042 (2020), eprint 2002.12534.
  * Athron et al. (2021) P. Athron, C. Balázs, D. H. Jacob, W. Kotlarski, D. Stöckinger, and H. Stöckinger-Kim (2021), eprint 2104.03691.
  * Baltz and Gondolo (2001) E. A. Baltz and P. Gondolo, Phys. Rev. Lett. 86, 5004 (2001), eprint hep-ph/0102147.
  * Everett et al. (2001) L. L. Everett, G. L. Kane, S. Rigolin, and L.-T. Wang, Phys. Rev. Lett. 86, 3484 (2001), eprint hep-ph/0102145.
  * Choudhury et al. (2001) D. Choudhury, B. Mukhopadhyaya, and S. Rakshit, Phys. Lett. B 507, 219 (2001), eprint hep-ph/0102199.
  * Cheung (2001) K.-m. Cheung, Phys. Rev. D 64, 033001 (2001), eprint hep-ph/0102238.
  * Calmet and Neronov (2002) X. Calmet and A. Neronov, Phys. Rev. D 65, 067702 (2002), eprint hep-ph/0104278.
  * Xiong and Yang (2001) Z.-H. Xiong and J. M. Yang, Phys. Lett. B 508, 295 (2001), eprint hep-ph/0102259.
  * Cox et al. (2021) P. Cox, C. Han, and T. T. Yanagida (2021), eprint 2104.03290.
  * Chakrabarty et al. (2018) N. Chakrabarty, C.-W. Chiang, T. Ohata, and K. Tsumura, JHEP 12, 104 (2018), eprint 1807.08167.
  * Chakrabarty (2020) N. Chakrabarty (2020), eprint 2010.05215.
  * Foot et al. (1994) R. Foot, X. G. He, H. Lew, and R. R. Volkas, Phys. Rev. D 50, 4571 (1994), eprint hep-ph/9401250.
  * Altmannshofer et al. (2016a) W. Altmannshofer, C.-Y. Chen, P. S. Bhupal Dev, and A. Soni, Phys. Lett. B 762, 389 (2016a), eprint 1607.06832.
  * Foldenauer and Jaeckel (2017) P. Foldenauer and J. Jaeckel, JHEP 05, 010 (2017), eprint 1612.07789.
  * Baek et al. (2001) S. Baek, N. G. Deshpande, X. G. He, and P. Ko, Phys. Rev. D 64, 055006 (2001), eprint hep-ph/0104141.
  * He et al. (1991a) X. G. He, G. C. Joshi, H. Lew, and R. R. Volkas, Phys. Rev. D 43, 22 (1991a).
  * He et al. (1991b) X.-G. He, G. C. Joshi, H. Lew, and R. R. Volkas, Phys. Rev. D 44, 2118 (1991b).
  * Bélanger et al. (2015) G. Bélanger, C. Delaunay, and S. Westhoff, Phys. Rev. D 92, 055021 (2015), eprint 1507.06660.
  * Heeck and Rodejohann (2011) J. Heeck and W. Rodejohann, Phys. Rev. D 84, 075007 (2011), eprint 1107.5238.
  * Lees et al. (2016) J. P. Lees et al. (BaBar), Phys. Rev. D 94, 011102 (2016), eprint 1606.03501.
  * Mishra et al. (1991) S. R. Mishra et al. (CCFR), Phys. Rev. Lett. 66, 3117 (1991).
  * Kaneta and Shimomura (2017) Y. Kaneta and T. Shimomura, PTEP 2017, 053B04 (2017), eprint 1701.00156.
  * Chun et al. (2019) E. J. Chun, A. Das, J. Kim, and J. Kim, JHEP 02, 093 (2019), [Erratum: JHEP 07, 024 (2019)], eprint 1811.04320.
  * Komatsu et al. (2011) E. Komatsu et al. (WMAP), Astrophys. J. Suppl. 192, 18 (2011), eprint 1001.4538.
  * Chiu (1966) H.-Y. Chiu, Phys. Rev. Lett. 17, 712 (1966).
  * Hall et al. (2010) L. J. Hall, K. Jedamzik, J. March-Russell, and S. M. West, JHEP 03, 080 (2010), eprint 0911.1120.
  * Bélanger et al. (2014) G. Bélanger, K. Kannike, A. Pukhov, and M. Raidal, JCAP 06, 021 (2014), eprint 1403.4960.
  * Baker et al. (2015) M. J. Baker et al., JHEP 12, 120 (2015), eprint 1510.03434.
  * Foldenauer (2019) P. Foldenauer, Phys. Rev. D 99, 035007 (2019), eprint 1808.03647.
  * Altmannshofer et al. (2016b) W. Altmannshofer, S. Gori, S. Profumo, and F. S. Queiroz, JHEP 12, 106 (2016b), eprint 1609.04026.
  * Biswas et al. (2017) A. Biswas, S. Choubey, and S. Khan, JHEP 02, 123 (2017), eprint 1612.03067.
  * Das et al. (2021) P. Das, M. K. Das, and N. Khan (2021), eprint 2104.03271.
  * Cheung et al. (2014) K. Cheung, Y.-L. S. Tsai, P.-Y. Tseng, T.-C. Yuan, and A. Zee, Nucl. Phys. B Proc. Suppl. 246-247, 116 (2014).
  * He et al. (2009) X.-G. He, T. Li, X.-Q. Li, J. Tandean, and H.-C. Tsai, Phys. Rev. D 79, 023521 (2009), eprint 0811.0658.
  * Chiang et al. (2013) C.-W. Chiang, T. Nomura, and J. Tandean, Phys. Rev. D 87, 073004 (2013), eprint 1205.6416.
  * Farzan and Akbarieh (2012) Y. Farzan and A. R. Akbarieh, JCAP 10, 026 (2012), eprint 1207.4272.
  * Farzan and Akbarieh (2014) Y. Farzan and A. R. Akbarieh, JCAP 11, 015 (2014), eprint 1408.2950.
  * Biswas et al. (2016a) A. Biswas, S. Choubey, and S. Khan, JHEP 09, 147 (2016a), eprint 1608.04194.
  * Chakrabarty et al. (2015) N. Chakrabarty, D. K. Ghosh, B. Mukhopadhyaya, and I. Saha, Phys. Rev. D 92, 015002 (2015), eprint 1501.03700.
  * Kannike (2016) K. Kannike, Eur. Phys. J. C 76, 324 (2016), [Erratum: Eur.Phys.J.C 78, 355 (2018)], eprint 1603.02680.
  * Biswas et al. (2016b) A. Biswas, S. Choubey, and S. Khan, JHEP 08, 114 (2016b), eprint 1604.06566.
  * Agostini et al. (2019) M. Agostini et al. (Borexino), Phys. Rev. D 100, 082004 (2019), eprint 1707.09279.
  * Abdullah et al. (2018) M. Abdullah, J. B. Dent, B. Dutta, G. L. Kane, S. Liao, and L. E. Strigari, Phys. Rev. D 98, 015005 (2018), eprint 1803.01224.
  * Altmannshofer et al. (2014) W. Altmannshofer, S. Gori, M. Pospelov, and I. Yavin, Phys. Rev. Lett. 113, 091801 (2014), eprint 1406.2332.
  * Tanabashi et al. (2018) M. Tanabashi et al. (Particle Data Group), Phys. Rev. D 98, 030001 (2018).
  * Kolb and Turner (1981) E. W. Kolb and M. S. Turner, Nature 294, 521 (1981).
  * Aghanim et al. (2020) N. Aghanim et al. (Planck), Astron. Astrophys. 641, A6 (2020), eprint 1807.06209.
  * Griest and Seckel (1991) K. Griest and D. Seckel, Phys. Rev. D 43, 3191 (1991).
  * Alloul et al. (2014) A. Alloul, N. D. Christensen, C. Degrande, C. Duhr, and B. Fuks, Comput. Phys. Commun. 185, 2250 (2014), eprint 1310.1921.
  * Bélanger et al. (2018) G. Bélanger, F. Boudjema, A. Goudelis, A. Pukhov, and B. Zaldivar, Comput. Phys. Commun. 231, 173 (2018), eprint 1801.03509.
  * Maurer (2016) V. Maurer, Comput. Phys. Commun. 198, 195 (2016), eprint 1503.01073.
  * Accomando et al. (2016) E. Accomando, C. Coriano, L. Delle Rose, J. Fiaschi, C. Marzo, and S. Moretti, JHEP 07, 086 (2016), eprint 1605.02910.
  * Robens and Stefaniak (2015) T. Robens and T. Stefaniak, Eur. Phys. J. C 75, 104 (2015), eprint 1501.02234.
  * López-Val and Robens (2014) D. López-Val and T. Robens, Phys. Rev. D 90, 114018 (2014), eprint 1406.1043.
  * Ahmed et al. (2010) Z. Ahmed et al. (CDMS-II), Science 327, 1619 (2010), eprint 0912.3592.
  * Aprile et al. (2017) E. Aprile et al. (XENON), Phys. Rev. Lett. 119, 181301 (2017), eprint 1705.06655.
  * Aprile et al. (2016) E. Aprile et al. (XENON), JCAP 04, 027 (2016), eprint 1512.07501.
  * Akerib et al. (2014) D. Akerib et al. (LUX), Phys. Rev. Lett. 112, 091303 (2014), eprint 1310.8214.
  * Abada et al. (2011) A. Abada, D. Ghaffor, and S. Nasri, Phys. Rev. D 83, 095021 (2011), eprint 1101.0365.
  * Basak et al. (2021) T. Basak, B. Coleppa, and K. Loho (2021), eprint 2105.09044.
  * Hoferichter et al. (2017) M. Hoferichter, P. Klos, J. Menéndez, and A. Schwenk, Phys. Rev. Lett. 119, 181803 (2017), eprint 1708.02245.