# Self-interacting Vectorial Dark Matter
in a SM-like Dark Sector

Van Que Tran  0000-0003-4643-4050<EMAIL_ADDRESS>Tsung-Dao Lee Institute &
School of Physics and Astronomy, Shanghai Jiao Tong University, Shanghai
200240, China  Faculty of Fundamental Sciences, PHENIKAA University, Yen
Nghia, Ha Dong, Hanoi 12116, Vietnam Thong T. Q. Nguyen  0000-0002-8460-0219
<EMAIL_ADDRESS>Institute of Physics, Academia Sinica, Nangang,
Taipei 11529, Taiwan Tzu-Chiang Yuan  0000-0001-8546-5031
<EMAIL_ADDRESS>Institute of Physics, Academia Sinica, Nangang,
Taipei 11529, Taiwan

###### Abstract

A $SU(2)_{D}\times U(1)_{D}$ gauge-Higgs sector, an exact dark copy of the
Standard Model (SM) one, is proposed. It is demonstrated that the dark gauge
bosons $\mathcal{W}^{(p,m)}$, in analogous to the SM $W^{\pm}$, can fulfill
the role as a self-interacting vector dark matter candidate, solving the core
versus cusp and missing satellites problems faced by the conventional paradigm
of collisionless weakly interacting massive particle. Constraints from
collider, astroparticle and cosmology on such a self-interacting vector dark
matter candidate are scrutinized. Implications for the future searches of
$\mathcal{W}^{(p,m)}$ in direct detection experiments are discussed.

##  I Introduction

The very conception of standard model of cosmology is the $\Lambda$CDM
paradigm in which the total mass-energy content of the Universe consists
approximately 26% of dark matter (DM), 5% of baryonic matter and the rest is
consistently described by a minuscule cosmological constant $\Lambda$. The DM
particles in $\Lambda$CDM cosmology is cold, collisionless and interact purely
gravitationally with each other and Standard Model (SM) particles. On large
scales greater than $\mathcal{O}({\rm Mpc})$, the structure of the Universe is
well described by the $\Lambda$CDM cosmology, as reflected in the predictions
of DM abundance, acoustic peaks in the microwave background radiation, matter
power spectrum, structure formation, etc. However, since 1990s various puzzles
had been posed on smaller scales for the cold dark matter (CDM) predictions.
Notably they are the core versus cusp, missing satellites and too-big-to-fail
problems. These are phrased as small scale “crisis” in the literature.

While SM does not have any DM candidate, many models beyond the SM (BSM) can
provide such a weakly-interacting massive particle (WIMP) to be the CDM. Most
popular models include the neutralino in the minimal supersymmetric standard
model (MSSM) motivated by the gauge hierarchy problem as well as the extremely
light axion particle associated with the solution of the strong CP problem in
QCD. Nevertheless all these DM candidates in the WIMP scenario suffer from the
small scale “crisis” mentioned in previous paragraph.

Over the years, many efforts were made to solve the small scale “crisis”
within the $\Lambda$CDM paradigm. There were (and still are) high hopes that
by including the baryon feedback processes such as gas cooling, star
formation, supernovae, and active galactic nuclei in the $N$-body simulations,
the “crisis” can be resolved. To the best of our knowledge, as of today, this
issue is not completely settled.

Alternatively one may want a paradigm shift by introducing self-interacting
dark matter (SIDM) first proposed by Spergel and Steinhardt Spergel:1999mh
some times ago as a solution to the core-cusp and missing satellites problems.
The idea is simple – SIDM is collisional and thus conducts heat. Self-
interactions can cause energy transfer from the hotter outer region to the
colder central region, thereby forming a core. As the core of a galaxy heats
up it expands, therefore decreasing the central density.

The simplest SIDM particle candidate is to augment the SM by adding a scalar
gauge singlet $\chi$ with a discrete symmetry $\chi\rightarrow-\chi$ that
allows only a quartic $\chi^{4}$ self-coupling and a quartic
$\chi^{2}\Phi_{H}^{\dagger}\Phi_{H}$ coupling where $\Phi_{H}$ is the SM Higgs
doublet as discussed in Silveira:1985rk ; Holz:2001cb ; Cheung:2012xb . At low
energy, the cross section $\sigma_{\chi\chi\rightarrow\chi\chi}$ is a constant
(velocity independent), dominated by the $\chi^{4}$ coupling while the $t$
channel exchange diagram due to the 125 GeV SM Higgs is negligible. To solve
the small scale “crisis” for dwarf galaxies, the SIDM cross section must be in
the following range Spergel:1999mh ,

$0.5\,{\rm cm^{2}/g}\lesssim\sigma/m_{\chi}\lesssim 6\,\mathrm{cm^{2}/g}\;,$
(1)

which implies Holz:2001cb

$m_{\chi}\sim\left(0.03-0.08\right)\lambda_{\chi}^{2/3}\,\mathrm{GeV}\;,$ (2)

where $\lambda_{\chi}$ is the $\chi^{4}$ coupling. Thus $\chi$ is a sub-GeV
SIDM if its self-coupling is of order one.

In recent years, there is a resurgence of interest in SIDM due to the series
of works in Tulin:2013teo ; Kaplinghat:2015aga ; Tulin:2017ara ;
Colquhoun:2020adl . Astrophysical data accumulated over the past several
decades for dwarf galaxies, low surface brightness (LSB) spiral galaxies and
galaxy clusters indicating DM self-interaction cross sections can have
significant velocity dispersion. The range in (1) has been extended to
Tulin:2017ara

$0.1\,{\rm cm^{2}/g}\lesssim\sigma/m_{\chi}\lesssim 10\,{\rm cm^{2}/g}\;,$ (3)

depending on the average velocity of the DM $\chi$ in the observed
astrophysical objects. We will return to this issue in later section. For
details on the small scale “crisis”, baryon feedback and SIDM, the readers are
redirected to a recent review by Tulin and Yu Tulin:2017ara .

Many BSM models can give rise to non-negligible $t$ channel exchange
contribution to SIDM cross section with mediator mass much lighter than the DM
mass, like the dark $U(1)$ force interacts with a fermionic DM in
Aboubrahim:2020lnr . Moreover the SIDM cross sections in this class of BSM
models have nontrivial velocity dependent of the DM and have sizes in the
desirable range given in (3).

Scalar Fields | $SU(2)_{L}$ | $U(1)_{Y}$ | $SU(2)_{D}$ | $U(1)_{D}$
---|---|---|---|---
$\Phi_{H}$ | $2$ | $\frac{1}{2}$ | $1$ | $0$
$\Phi_{D}$ | $1$ | $0$ | $2$ | $\frac{1}{2}$
Table 1: Scalar field contents in the model. Our convention is $Q=T^{3}+Y$
where $Q$ is the electric charge in unit of $e>0$, $T^{3}$ is the third
component of the generators of $SU(2)_{L}$ and $Y$ is the hypercharge.

In this paper, motivated by our earlier works on gauged two-Higgs-doublet
model (G2HDM) Huang:2015wts  and its minimal version Ramos:2021txu ;
Ramos:2021omo , we consider a $SU(2)_{D}\times U(1)_{D}$ gauge-Higgs sector
which is exactly a dark copy of the SM one $SU(2)_{L}\times U(1)_{Y}$. The
dark gauge bosons $\mathcal{W}^{(p,m)}$ of $SU(2)_{D}$, in parallel with the
SM $W^{\pm}$, can be a vectorial dark matter candidate. The main motivation or
advantage is that, due to the non-abelian cubic and quartic couplings,
$\mathcal{W}^{(p,m)}$ is naturally self-interacting. Furthermore, since
$\mathcal{W}^{(p,m)}$ is complex, one can have either symmetric (with both
$\mathcal{W}^{p}$ and $\mathcal{W}^{m}$) or asymmetric (with either
$\mathcal{W}^{p}$ or $\mathcal{W}^{m}$) DM scenario Graesser:2011wi ;
Iminniyaz:2011yp ; Lin:2011gj ; Zurek:2013wia . We will focus on the symmetric
case in what follows. The fermion content is the same as in SM, while the
scalar fields used in the model are shown in Table 1. The communication
between the SM and dark sectors are via three portals:

  * •

Kinetic mixing of the two abelian groups, the hypercharge $U(1)_{Y}$ in SM and
$U(1)_{D}$ in the dark sector,

  * •

Higgs mixing among the SM Higgs doublet $\Phi_{H}$ and its dark partner
$\Phi_{D}$, and

  * •

Mass mixings among the three neutral gauge bosons.

A Stueckelberg mass is introduced for $U(1)_{D}$ such that all gauge bosons
are massive except the SM photon. We will be interested SIDM in the sub-GeV to
TeV mass range with light mediators in order to resolve some of the well-known
issues like core-cusp and too-big-to-fail problems Tulin:2017ara  that one
faces in the more conventional WIMP CDM scenario. Indeed, the model consists
of one dark photon $\gamma^{\prime}$, one dark $Z^{\prime}$ and one dark
Higgs, all of which can play the role of a mediator $X$ for SIDM. These extra
particles can have important impacts in collider physics as well.

Computation of the SIDM cross section is highly nontrivial as pointed out in
Colquhoun:2020adl . Methods used for the computation depend on two
dimensionless parameters $\kappa$ and $\beta$ constructed from the following 4
quantities: coupling strength ($\alpha_{X}$) between the SIDM ($\mathcal{W}$)
and its mediator ($X$), their masses ($m_{\cal W}$ and $m_{X}$) and the
velocity $v$ of the DM,

$\kappa=\frac{m_{\cal
W}v}{m_{X}}\;,\quad\quad\beta=\frac{2\alpha_{X}m_{X}}{m_{\cal W}v^{2}}\;.$ (4)

The naïve Born approximation ($2\beta\kappa^{2}\ll 1$) based on perturbative
expansion in $\alpha_{X}$ would break down when the mediator mass is much
lighter than the DM mass resulting a long range force. Non-perturbative
effects must be taken into account in the semi-classical ($\kappa\geq 1$) and
quantum regimes ($\kappa\ll 1$), as suggested in Colquhoun:2020adl  which we
will follow closely in this work.

The content of this paper is unfolded as follows. In section II, we present
the Lagrangian of the gauge-Higgs sector in the model and discuss its mass
spectra after spontaneously symmetry breaking. Phenomenological constraints
for the vectorial DM $\mathcal{W}^{(p,m)}$ from collider physics, electroweak
precision tests, dark photon and dark $Z^{\prime}$ physics, relic density and
direct detection are considered in section III. Implementation of
$\mathcal{W}^{(p,m)}$ as SIDM is discussed in section IV, followed by our
numerical analysis in section V. We conclude in section VI. Details of the
relevant couplings in the neutral current interaction Lagrangian for the
fermions and gauge bosons are given in appendix A. Results of an alternative
mass hierarchy of the dark photon and dark $Z^{\prime}$ is given in appendix
B.

##  II The $SU(2)_{L}\times U(1)_{Y}\times SU(2)_{D}\times U(1)_{D}$ Model

As alluded to in previous section, the main advantage of having a non-abelian
gauge-Higgs dark sector of $SU(2)_{D}\times U(1)_{D}$ is that the extra gauge
bosons $\mathcal{W}^{(p,m)}$ can be naturally a self-interacting dark matter
candidate without the need of introducing other matter field. Thus the
Lagrangian of the gauge sector is given by

$\displaystyle{\mathcal{L}_{G}}={}$
$\displaystyle-\frac{1}{4}F^{a}_{\mu\nu}F^{a{}\mu\nu}-\frac{1}{4}B_{\mu\nu}B^{\mu\nu}-\frac{1}{4}\mathcal{F}^{a}_{\mu\nu}\mathcal{F}^{a}_{\mu\nu}-\frac{1}{4}\mathcal{X}_{\mu\nu}\mathcal{X}^{\mu\nu}$
$\displaystyle-\frac{1}{2}\epsilon\mathcal{X}_{\mu\nu}B^{\mu\nu}+\frac{1}{2}\left(\partial_{\mu}\sigma+M_{\mathcal{X}}{\mathcal{X}}_{\mu}\right)^{2}+\mathrm{gauge\;fixings}\;,$
(5)

where $F^{a}_{\mu\nu}$ and $\mathcal{F}^{a}_{\mu\nu}$ $(a=1,2,3)$ are the
$SU(2)_{L}$ and $SU(2)_{D}$ field strength tensors respectively,
$M_{\mathcal{X}}$ is the Stueckelberg mass Kors:2004dx ; Kors:2005uz ;
Cheung:2007ut  for the $U(1)_{D}$ gauge field ${\mathcal{X}}_{\mu}$ associated
with its field strength $\mathcal{X}_{\mu\nu}$, and $\epsilon$ is the kinetic
mixing parameter Holdom:1985ag ; Feldman:2007wj  between the SM hypercharge
field strength $B_{\mu\nu}$ and $\mathcal{X}_{\mu\nu}$. $\sigma$ is the
auxiliary scalar field introduced to implement the Stueckelberg mechanism to
provide a mass for the $U(1)_{D}$ gauge boson in a gauge invariant fashion.
The SM fermions are the same and they are inert under the dark group
$SU(2)_{D}\times U(1)_{D}$. Dark fermions are also possible but we do not
introduce them to keep thing as minimal as possible.

On the other hand, the Lagrangian for the scalar sector is given by

$\displaystyle{\mathcal{L}}_{S}={}$
$\displaystyle|D_{\mu}\Phi_{H}|^{2}+|D_{\mu}\Phi_{D}|^{2}-V\left(\Phi_{H},\Phi_{D}\right)\;,$
(6)

where the covariant derivatives of $D_{\mu}\Phi$ and $D_{\mu}\Phi_{D}$ are
given by the familiar expressions 111Here for $\Phi_{H}$ we have
$T^{\pm}=T^{1}\pm iT^{2}=\left(\sigma^{1}\pm i\sigma^{2}\right)/2$ and
$T^{3}=\sigma^{3}/2$ where $\sigma^{i}(i=1,2,3)$ are the Pauli matrices.
Similarly for $\Phi_{D}$ we have $\mathcal{T}^{\pm}=\mathcal{T}^{1}\pm
i\mathcal{T}^{2}=\left(\sigma^{1}\pm i\sigma^{2}\right)/2$ and
$\mathcal{T}^{3}=\sigma^{3}/2$. $W^{\pm}_{\mu}=\left(W^{1}_{\mu}\mp
iW^{2}_{\mu}\right)/\sqrt{2}$ and $W^{3}_{\mu}$ are the gauge fields of
$SU(2)_{L}$, $\mathcal{W}^{(p,m)}_{\mu}=\left(\mathcal{W}^{1}_{\mu}\mp
i\mathcal{W}^{2}_{\mu}\right)/\sqrt{2}$ and $\mathcal{W}^{3}_{\mu}$ are the
gauge fields of $SU(2)_{D}$, and finally $B_{\mu}$ and $\mathcal{X}_{\mu}$ are
the gauge fields of $U(1)_{Y}$ and $U(1)_{D}$ respectively.

$\displaystyle D_{\mu}\Phi_{H}=$
$\displaystyle\left[\partial_{\mu}-i\frac{g}{2}\left(W^{+}_{\mu}T^{+}+W^{-}_{\mu}T^{-}\right)-igW_{\mu}^{3}T^{3}-ig^{\prime}\frac{1}{2}B_{\mu}\right]\Phi_{H}\;,$
$\displaystyle D_{\mu}\Phi_{D}=$
$\displaystyle\left[\partial_{\mu}-i\frac{g_{D}}{2}\left({\mathcal{W}}^{p}_{\mu}{\mathcal{T}}^{p}+{\mathcal{W}}^{m}_{\mu}{\mathcal{T}}^{m}\right)-ig_{D}{\mathcal{W}}^{3}_{\mu}\mathcal{T}^{3}-ig_{D}^{\prime}\frac{1}{2}{\mathcal{X}}_{\mu}\right]\Phi_{D}\;;$
(7)

and $V(\Phi_{H},\Phi_{D})$ is the most general Higgs potential invariant under
both $SU(2)_{L}\times U(1)_{Y}$ and $SU(2)_{D}\times U(1)_{D}$

$\displaystyle V\left(\Phi_{H},\Phi_{D}\right)={}$
$\displaystyle-\mu^{2}\Phi_{H}^{\dagger}\Phi_{H}+\lambda\left(\Phi_{H}^{\dagger}\Phi_{H}\right)^{2}-\mu^{2}_{D}\Phi_{D}^{\dagger}\Phi_{D}+\lambda_{D}\left(\Phi_{D}^{\dagger}\Phi_{D}\right)^{2}$
$\displaystyle+\lambda^{\prime}\left(\Phi_{H}^{\dagger}\Phi_{H}\right)\left(\Phi_{D}^{\dagger}\Phi_{D}\right)\;.$
(8)

Neutral current interaction Lagrangian for the fermions and gauge bosons can
be found in appendix A.

###  II.1 Spontaneous Symmetry Breaking

To study the spontaneous symmetry breaking (SSB) in the model, we decompose
the fields in the two scalar doublets according to

$\displaystyle\Phi_{H}$ $\displaystyle=$ $\displaystyle\begin{pmatrix}G^{+}\\\
\frac{1}{\sqrt{2}}\left(v_{H}+h+iG^{0}\right)\end{pmatrix}\;,\;\;\Phi_{D}=\begin{pmatrix}G^{p}\\\
\frac{1}{\sqrt{2}}\left(v_{D}+h_{D}+iG^{0}_{D}\right)\end{pmatrix}\;.$ (9)

Here $v_{H}=246$ GeV is the SM vacuum expectation value (VEV), and $v_{D}$ is
the unknown hidden VEV,

$\displaystyle\langle\Phi_{H}\rangle=\begin{pmatrix}0\\\
\frac{v_{H}}{\sqrt{2}}\end{pmatrix},\;\langle\Phi_{D}\rangle=\begin{pmatrix}0\\\
\frac{v_{D}}{\sqrt{2}}\end{pmatrix}.\;$ (10)

The potential for the VEVs is then

$V(v_{H},v_{D})=\frac{1}{4}\left(-2\mu^{2}v_{H}^{2}+\lambda
v_{H}^{4}-2\mu_{D}^{2}v_{D}^{2}+\lambda_{D}v_{D}^{4}+\lambda^{\prime}v_{H}^{2}v_{D}^{2}\right)\;.$
(11)

To determine the true vacuum, we impose the tadpole conditions $\partial
V/\partial v_{H}=0$ and $\partial V/\partial v_{D}=0$, which lead to

$\displaystyle v_{H}\left(-2\mu^{2}+2\lambda
v_{H}^{2}+\lambda^{\prime}v_{D}^{2}\right)=0\;,$ (12) $\displaystyle
v_{D}\left(-2\mu_{D}^{2}+2\lambda_{D}v_{D}^{2}+\lambda^{\prime}v_{H}^{2}\right)=0\;.$
(13)

Non-trivial solutions for the VEVs are

$\displaystyle v_{H}^{2}=$
$\displaystyle\frac{4\lambda_{D}\mu^{2}-2\lambda^{\prime}\mu_{D}^{2}}{4\lambda\lambda_{D}-\lambda^{\prime
2}}\;,$ (14) $\displaystyle v_{D}^{2}=$
$\displaystyle\frac{4\lambda\mu_{D}^{2}-2\lambda^{\prime}\mu^{2}}{4\lambda\lambda_{D}-\lambda^{\prime
2}}\;.$ (15)

Besides requiring the vacuum to be non-tachyonic $v_{H}^{2}>0$ and
$v_{D}^{2}>0$, to trust our perturbative calculations, we must impose

$|\lambda|,\;|\lambda_{D}|\;\textrm{and}\;|\lambda^{\prime}|\leq 8\pi\;.$ (16)

Tree level perturbative unitarity in the scalar sector implies further that
Arhrib:2012ia ; Arhrib:2013ela

$|\lambda|\;\textrm{and}\;|\lambda_{D}|\leq\frac{4\pi}{3}\;.$ (17)

Moreover, requiring the potential bounded from below implies Branco:2011iw

$\lambda>0,\;\lambda_{D}>0\;\textrm{and}\;2\sqrt{\lambda\lambda_{D}}+\lambda^{\prime}>0\;.$
(18)

###  II.2 Mass Spectra

####  II.2.1 Scalar Bosons

We will work in the ’t Hooft-Landau gauge in which the Goldstone bosons from
both SM and dark sectors are massless (in view of the above tadpole conditions
(12) and (13))

$\displaystyle m_{G^{\pm}}^{2}$ $\displaystyle=$ $\displaystyle
m_{G^{0}}^{2}=-\mu^{2}+\lambda
v_{H}^{2}+\frac{1}{2}\lambda^{\prime}v_{D}^{2}=0\;,$ (19) $\displaystyle
m_{G^{(p,m)}}^{2}$ $\displaystyle=$ $\displaystyle
m_{G^{0}_{D}}^{2}=-\mu_{D}^{2}+\lambda_{D}v_{D}^{2}+\frac{1}{2}\lambda^{\prime}v_{H}^{2}=0\;.$
(20)

The remaining two components $h$ and $h_{D}$ from the two doublets $\Phi_{H}$
and $\Phi_{D}$ mix according to

${\mathcal{M}}^{2}_{S}=\begin{pmatrix}-\mu^{2}+3\lambda
v_{H}^{2}+\frac{1}{2}\lambda^{\prime}v_{D}^{2}&\lambda^{\prime}v_{H}v_{D}\\\
\lambda^{\prime}v_{H}v_{D}&-\mu_{D}^{2}+3\lambda_{D}v_{D}^{2}+\frac{1}{2}\lambda^{\prime}v_{H}^{2}\end{pmatrix}=\begin{pmatrix}2\lambda
v_{H}^{2}&\lambda^{\prime}v_{H}v_{D}\\\
\lambda^{\prime}v_{H}v_{D}&2\lambda_{D}v_{D}^{2}\end{pmatrix}\;.$ (21)

In the last equality we have used again the tadpole conditions to eliminate
$\mu^{2}$ and $\mu_{D}^{2}$. The physical states $h_{1}$ and $h_{2}$ are
related to $h$ and $h_{D}$ by

$\displaystyle h=$ $\displaystyle h_{1}\cos\alpha-h_{2}\sin\alpha\;,$ (22)
$\displaystyle h_{D}=$ $\displaystyle h_{1}\sin\alpha+h_{2}\cos\alpha\;,$ (23)

with the mixing angle

$\sin 2\alpha=\frac{\lambda^{\prime}v_{H}v_{D}}{\left[\left(\lambda
v_{H}^{2}-\lambda_{D}v_{D}^{2}\right)^{2}+\lambda^{\prime
2}v_{H}^{2}v_{D}^{2}\right]^{1/2}}\;,$ (24)

and their masses are given by

$m_{h_{1},h_{2}}^{2}=\left(\lambda
v_{H}^{2}+\lambda_{D}v_{D}^{2}\right)\mp\left[\left(\lambda
v_{H}^{2}-\lambda_{D}v_{D}^{2}\right)^{2}+\lambda^{\prime
2}v_{H}^{2}v_{D}^{2}\right]^{1/2}\;.$ (25)

We will identify the lighter Higgs $h_{1}$ with mass $m_{h_{1}}=125.25\pm
0.17$ GeV PDG2022  corresponds to the scalar Higgs boson observed at LHC.

One can invert the above formulas (24) and (25) to trade the fundamental
parameters $\lambda$, $\lambda_{D}$ and $\lambda^{\prime}$ with the two
physical masses $m_{1}^{2}$, $m_{2}^{2}$ and the mixing angle $\alpha$, namely

$\displaystyle\lambda=$
$\displaystyle\frac{1}{2v_{H}^{2}}\left(m_{h_{1}}^{2}\cos^{2}\alpha+m_{h_{2}}^{2}\sin^{2}\alpha\right)\;,$
(26) $\displaystyle\lambda_{D}=$
$\displaystyle\frac{1}{2v_{D}^{2}}\left(m_{h_{1}}^{2}\sin^{2}\alpha+m_{h_{2}}^{2}\cos^{2}\alpha\right)\;,$
(27) $\displaystyle\lambda^{\prime}=$
$\displaystyle\frac{1}{2v_{H}v_{D}}\left[\left(m_{h_{2}}^{2}-m_{h_{1}}^{2}\right)\sin
2\alpha\right]\;.$ (28)

These equations are useful for numerical studies.

####  II.2.2 Gauge Bosons

The mass spectrum of gauge bosons from spontaneous symmetry breaking can be
obtained by looking at the terms in $|D_{\mu}\Phi_{H}|^{2}$ and
$|D_{\mu}\Phi_{D}|^{2}$ which are second order in the VEVs ($v_{H},v_{D}$) as
well as second order in the gauge couplings
($g,g^{\prime},g_{D},g_{D}^{\prime}$).

There is only one charged vector boson, the SM $W^{\pm}$, with mass given by
the standard formula

$m_{W^{\pm}}=\frac{1}{2}gv_{H}\;.$ (29)

Experimentally we have $m_{W}=80.377\pm 0.012\;{\rm GeV}$ PDG2022 .

The complex dark gauge boson ${\mathcal{W}}^{p}$
(${\mathcal{W}}^{m}=({\mathcal{W}}^{p})^{*}$) carries one unit of $+(-)$ dark
charge but no electric charge. In analogous to (29), their masses are simply

$m_{\mathcal{W}^{(p,m)}}=\frac{1}{2}g_{D}v_{D}\;.$ (30)

All the other 4 gauge bosons $W^{3},B,{\mathcal{W}^{3}}$ and $\mathcal{X}$ are
electrically neutral and in general mixed together. In the basis of

$V_{N}=\begin{pmatrix}B\\\ W^{3}\\\ \mathcal{X}\\\
{\mathcal{W}}^{3}\end{pmatrix}\;,$ (31)

one finds

$\displaystyle{\mathcal{M}}_{N}^{2}$ $\displaystyle=$
$\displaystyle\begin{pmatrix}\frac{1}{4}g^{\prime\,2}v_{H}^{2}&-\frac{1}{4}gg^{\prime}v_{H}^{2}&0&0\\\
-\frac{1}{4}gg^{\prime}v_{H}^{2}&\frac{1}{4}g^{2}v_{H}^{2}&0&0\\\
0&0&\frac{1}{4}g_{D}^{\prime
2}v_{D}^{2}+M_{\mathcal{X}}^{2}&-\frac{1}{4}g_{D}g_{D}^{\prime}v_{D}^{2}\\\
0&0&-\frac{1}{4}g_{D}g_{D}^{\prime}v_{D}^{2}&\frac{1}{4}g_{D}^{2}v_{D}^{2}\end{pmatrix}\;.$
(32)

Note that we have included a Stueckelberg mass $M_{\mathcal{X}}$ for the
$U(1)_{D}$, while for the $U(1)_{Y}$ we set its Stueckelberg mass $M_{Y}$ to
be zero.

Before diagonalize the above mass matrix ${\mathcal{M}}_{N}^{2}$, we need to
first diagonalize the kinetic terms for $B_{\mu\nu}$ and
$\mathcal{X}_{\mu\nu}$ due to their kinetic mixing. This can be accomplished
by performing a $GL(4)$ transformation $K$ given by Feldman:2007wj

$K=\begin{pmatrix}1&0&-S_{\epsilon}&0\\\ 0&1&0&0\\\ 0&0&C_{\epsilon}&0\\\
0&0&0&1\end{pmatrix}\;,\;\;\;{\rm
with}\;C_{\epsilon}=\frac{1}{\sqrt{1-\epsilon^{2}}}\;{\rm
and}\;S_{\epsilon}=\epsilon C_{\epsilon}\;.$ (33)

This induces a change in the mass matrix
${\mathcal{M}}_{N}^{2}\rightarrow{\mathcal{M}}_{N}^{\prime
2}=K^{T}\cdot{\mathcal{M}}_{N}^{2}\cdot K$ and $V_{N}\rightarrow
V^{\prime}_{N}=K^{-1}\cdot V_{N}$, i.e.

$\displaystyle{\mathcal{M}}_{N}^{\prime
2}=\begin{pmatrix}\frac{1}{4}g^{\prime\,2}v_{H}^{2}&-\frac{1}{4}gg^{\prime}v_{H}^{2}&-\frac{1}{4}g^{\prime
2}v_{H}^{2}S_{\epsilon}&0\\\
-\frac{1}{4}gg^{\prime}v_{H}^{2}&\frac{1}{4}g^{2}v_{H}^{2}&\frac{1}{4}gg^{\prime}v_{H}^{2}S_{\epsilon}&0\\\
-\frac{1}{4}g^{\prime
2}v_{H}^{2}S_{\epsilon}&\frac{1}{4}gg^{\prime}v_{H}^{2}S_{\epsilon}&\left(\frac{1}{4}g_{D}^{\prime
2}v_{D}^{2}+M_{\mathcal{X}}^{2}\right)C^{2}_{\epsilon}+\frac{1}{4}g^{\prime
2}v_{H}^{2}S^{2}_{\epsilon}&-\frac{1}{4}g_{D}g_{D}^{\prime}v_{D}^{2}C_{\epsilon}\\\
0&0&-\frac{1}{4}g_{D}g_{D}^{\prime}v_{D}^{2}C_{\epsilon}&\frac{1}{4}g_{D}^{2}v_{D}^{2}\end{pmatrix}\;,$
(34)

and

$V^{\prime}_{N}\equiv\begin{pmatrix}B^{\prime}\\\ W^{3\prime}\\\
\mathcal{X}^{\prime}\\\
{\mathcal{W}}^{3\prime}\end{pmatrix}=\begin{pmatrix}B+\epsilon\mathcal{X}\\\
W^{3}\\\ \frac{1}{C_{\epsilon}}\mathcal{X}\\\
{\mathcal{W}}^{3}\end{pmatrix}\;.$ (35)

Despite its unappealing look, (34) has zero determinant so it has at least one
zero eigenvalue 222If $M_{\mathcal{X}}$ vanishes, (34) will have two zero
eigenvalues. Besides the massless photon in the SM, one would have a massless
dark photon as well.. Indeed, by using the following orthogonal transformation
$V_{N}^{\prime}\rightarrow V_{N}^{\prime\prime}=O_{W}^{T}\cdot V^{\prime}$
where

$O_{W}=\begin{pmatrix}c_{W}&-s_{W}&0&0\\\ s_{W}&c_{W}&0&0\\\ 0&0&1&0\\\
0&0&0&1\end{pmatrix}\;\;\;{\rm
with}\;\;\left(c_{W},s_{W}\right)=\frac{1}{\sqrt{g^{2}+g^{\prime\,2}}}\left(g,g^{\prime}\right)\;,$
(36)

and hence 333Here $A_{\rm SM}=c_{W}B+s_{W}W^{3}$ and $Z_{\rm
SM}=-s_{W}B+c_{W}W^{3}$ are the SM expressions.

$V^{\prime\prime}_{N}\equiv\begin{pmatrix}A\\\ Z\\\
\mathcal{X}^{\prime\prime}\\\
{\mathcal{W}}^{3\prime\prime}\end{pmatrix}=\begin{pmatrix}A_{\rm
SM}+c_{W}\epsilon\mathcal{X}\\\ Z_{\rm SM}-s_{W}\epsilon\mathcal{X}\\\
\frac{1}{C_{\epsilon}}\mathcal{X}\\\ {\mathcal{W}}^{3}\end{pmatrix}\;,$ (37)

one can bring ${\mathcal{M}}_{N}^{\prime\,2}$ into the following form

$O_{W}^{T}\cdot{\mathcal{M}}_{N}^{\prime\,2}\cdot O_{W}=\begin{pmatrix}0&{\bf
0}^{\rm T}\\\ {\bf 0}&\hat{\mathcal{M}}_{N}^{2}\end{pmatrix}\;,$ (38)

where $\hat{\mathcal{M}}_{N}^{2}$ is a reduced $3\times 3$ matrix given by

$\hat{\mathcal{M}}_{N}^{2}=\begin{pmatrix}\frac{1}{4}\left(g^{2}+g^{\prime
2}\right)v_{H}^{2}&\frac{1}{4}g^{\prime}\sqrt{g^{2}+g^{\prime
2}}v_{H}^{2}S_{\epsilon}&0\\\ \frac{1}{4}g^{\prime}\sqrt{g^{2}+g^{\prime
2}}v_{H}^{2}S_{\epsilon}&\left(\frac{1}{4}g_{D}^{\prime
2}v_{D}^{2}+M_{\mathcal{X}}^{2}\right)C^{2}_{\epsilon}+\frac{1}{4}g^{\prime
2}v_{H}^{2}S_{\epsilon}^{2}&-\frac{1}{4}g_{D}g_{D}^{\prime}v_{D}^{2}C_{\epsilon}\\\
0&-\frac{1}{4}g_{D}g_{D}^{\prime}v_{D}^{2}C_{\epsilon}&\frac{1}{4}g_{D}^{2}v_{D}^{2}\end{pmatrix}\;.$
(39)

One can see from (37) and (38) that the $A$ field has zero mass and can be
identified as the physical photon in the model, while the $Z$ field will be
further mixed with $\mathcal{X}^{\prime\prime}$ and
${\mathcal{W}}^{3\prime\prime}$ with the mass matrix defined in (39). We thus
perform a final orthogonal transformation $\mathcal{O}$ to diagonalize
$\hat{\mathcal{M}}_{N}^{2}$ as follows

$\begin{pmatrix}Z\\\ \mathcal{X}^{\prime\prime}\\\
\mathcal{W}^{3\prime\prime}\end{pmatrix}\rightarrow\begin{pmatrix}Z_{1}\\\
Z_{2}\\\ Z_{3}\end{pmatrix}=\mathcal{O}^{T}\cdot\begin{pmatrix}Z\\\
\mathcal{X}^{\prime\prime}\\\ \mathcal{W}^{3\prime\prime}\end{pmatrix}\;,$
(40)

and

$\hat{\mathcal{M}}_{N}^{2}\rightarrow\mathcal{O}^{T}\cdot\hat{\mathcal{M}}_{N}^{2}\cdot\mathcal{O}=\begin{pmatrix}m_{1}^{2}&0&0\\\
0&m_{2}^{2}&0\\\ 0&0&m_{3}^{2}\end{pmatrix}\;,$ (41)

where $m_{i}^{2}$ is the mass-squared of the physical field $Z_{i}$ for
$i=1,2,3$. Note that ${\rm
Det}\hat{\mathcal{M}}_{N}^{2}=g_{D}^{2}(g^{2}+g_{D}^{2})v_{H}^{2}v_{D}^{2}C_{\epsilon}^{2}M_{\mathcal{X}}^{2}/16$.
In the limit of $\epsilon\rightarrow 0$ and $M_{\mathcal{X}}\rightarrow 0$, it
is clear from the structure of $\hat{\mathcal{M}}_{N}^{2}$ in (39) that
$Z_{1}\rightarrow Z_{\rm SM}$ with $m_{1}^{2}\rightarrow m^{2}_{Z_{\rm
SM}}=(g^{2}+g^{\prime 2})^{2}v_{H}^{2}/4$; while
$Z_{2}\rightarrow\gamma_{D}=c_{D}\mathcal{\chi}^{\prime\prime}+s_{D}\mathcal{W}^{3\prime\prime}$
and
$Z_{3}\rightarrow\mathcal{Z}_{D}=-s_{D}\mathcal{\chi}^{\prime\prime}+c_{D}\mathcal{W}^{3\prime\prime}$
with $m^{2}_{2}\rightarrow m^{2}_{\gamma_{D}}=0$ and $m^{2}_{3}\rightarrow
m^{2}_{\mathcal{Z}_{D}}=m^{2}_{\mathcal{W}^{(p,m)}}/c^{2}_{D}$, where
$(c_{D},s_{D})=(g_{D},g_{D}^{\prime})/\sqrt{g_{D}^{2}+g_{D}^{\prime 2}}$ in
analogous to the SM Weinberg angle for the dark sector. In order to have a
successful symmetry breaking pattern for general values of $\epsilon$ and
$M_{\mathcal{X}}$, it is necessary for $\hat{\mathcal{M}}_{N}^{2}$ to provide
an eigenvalue with a square root within the experimental range of the physical
$Z$ mass $m_{Z}=91.1876\pm 0.0021$ GeV PDG2022 . While it is possible to
obtain analytic formulas for the rotation angles in $\mathcal{O}$ in terms of
the fundamental parameters in the Lagrangian as was done in previous work
Huang:2019obt , we will diagonalize $\hat{\mathcal{M}}_{N}^{2}$ numerically in
this work.

##  III Phenomenological Constraints

As mentioned in the Introduction section, the communication between the SM and
dark sectors is controlled by the kinetic mixing parameter $\epsilon$, the
Higgs mixing parameter $\lambda^{\prime}$, and the mass mixings among the
three neutral gauge bosons. By tuning these mixing effects sufficiently small,
the SM physics remains intact. The other new and relevant parameters in the
model are $g_{D}$, $g_{D}^{\prime}$, $v_{D}$ and $\lambda_{D}$, which we would
like to constrain by collider physics, astrophysics and cosmology in this
work.

###  III.1 Collider Physics

The extension of the scalar sector in the model can be constrained by the data
from Higgs search experiments and the measurements of the SM-like Higgs boson
at LEP, Tevatron and LHC.

Here, the Higgs phenomenology at colliders is similar to that of one real
singlet scalar extension of the SM. The most relevant parameters are the
mixing angle $\alpha$ and the mass of the extra Higgs boson $m_{h_{2}}$. The
couplings of the physical scalars $h_{1}$ and $h_{2}$ to the SM particles can
be given as

${\cal L}_{\rm Higgs}\supset\frac{h_{1}\cos\alpha-
h_{2}\sin\alpha}{v_{H}}\left(2m_{W}^{2}W^{+}_{\mu}W^{-\mu}+m_{Z}Z_{\mu}Z^{\mu}-\sum_{f}m_{f}{\bar{f}}f\right)\;.$
(42)

One can see that the SM-like Higgs boson coupling to the SM particles are
modified by a factor of $\cos\alpha$. The Higgs boson signal strength then can
be given as

$\mu_{h_{1}}\equiv\cos^{2}\alpha\frac{{\rm BR}(h_{1}\rightarrow{\rm SM})}{{\rm
BR^{SM}}(h_{1}\rightarrow{\rm SM})}\;,$ (43)

where $\text{BR}^{\text{SM}}(h_{1}\rightarrow\text{SM})\equiv 1$ and
$\text{BR}(h_{1}\rightarrow\text{SM})=\frac{\Gamma_{h_{1}}^{\text{SM}}\cos^{2}\alpha}{\Gamma_{h_{1}}^{\text{SM}}\cos^{2}\alpha+\Gamma_{h_{1}}^{\text{DS}}}$
with $\Gamma_{h_{1}}^{\text{DS}}$ is the partial decay width of $h_{1}$ to
dark sector particles. Therefore, one obtains

$\mu_{h_{1}}=\frac{\Gamma_{h_{1}}^{\text{SM}}\cos^{4}\alpha}{\Gamma_{h_{1}}^{\text{SM}}\cos^{2}\alpha+\Gamma_{h_{1}}^{\text{DS}}}\;.$
(44)

It is expected that $\Gamma_{h_{1}}^{\text{DS}}\ll\Gamma_{h_{1}}^{\text{SM}}$
thus the signal strength simply scales as $\cos^{2}\alpha$. Using the current
Higgs signal strengths measurement $\mu_{h}=1.05\pm 0.06$ by ATLAS
ATLAS:2022vkf , one can obtain a bound on the mixing angle
$|\sin\alpha|\lesssim 0.2$ at $95\%$ C.L.

If $m_{h_{1}}\geq 2m_{\cal W}$, the Higgs boson $h_{1}$ can decay invisibly
into a pair of ${\cal W}$. The invisible branching ratio can be given by

${\rm BR}({h_{1}\rightarrow{\rm inv}})=\frac{\Gamma({h_{1}\rightarrow{\cal
W}^{p}{\cal
W}^{m}})}{\Gamma_{h_{1}}^{\text{SM}}\cos^{2}\alpha+\Gamma_{h_{1}}^{\text{DS}}}\;,$
(45)

where the invisible decay width can be given as

$\Gamma({h_{1}\rightarrow{\cal W}^{p}{\cal
W}^{m}})=\frac{g_{D}^{4}\left(v_{D}\sin\alpha\right)^{2}}{256\pi}\frac{m_{h_{1}}^{3}}{m_{{\cal
W}}^{4}}\left(1-\tau_{{\cal W}}+\frac{3}{4}\tau_{{\cal
W}}^{2}\right)\sqrt{1-\tau_{{\cal W}}}\;,$ (46)

with $\tau_{{\cal W}}=4m_{{\cal W}}^{2}/m_{h_{1}}^{2}$. Recently, ATLAS sets a
limit of ${\rm BR}({h_{1}\rightarrow{\rm inv}})<0.13$ at $95\%$ C.L., assuming
that the Higgs boson production cross section via vector boson fusion is
comparable to the SM prediction ATLAS:2022vkf .

Figure 1: Upper bounds on the mixing angle $\alpha$ as a function of the heavy
Higgs mass $m_{h_{2}}$. The color shaded regions represent the exclusion
regions from the Higgs signal strength measurements at ATLAS ATLAS:2022vkf
(gray region), the oblique parameter constraint (blue region) and the di-boson
searches via gluon-gluon fusion (purple region) and vector boson fusion (green
region) channels at ATLAS ATLAS:2020fry .

The heavy Higgs boson $h_{2}$ can interact with the SM particles through the
mixing with the SM-like Higgs boson $h_{1}$. If kinematically allowed, $h_{2}$
can decay to a pair of massive gauge bosons, a pair of fermions, a pair of SM-
like Higgs bosons or a pair of dark sector particles. The resonance production
rate of $h_{2}$ to a pair of SM particles $X$ can be given as

$\sigma(pp\rightarrow XX)\equiv\sigma(pp\rightarrow h_{2})\times{\rm
BR}(h_{2}\rightarrow XX)\;,$ (47)

where $\sigma(pp\rightarrow h_{2})=\sin^{2}\alpha\,\sigma^{\rm
SM}(pp\rightarrow h)|_{m_{h}=m_{h_{2}}}$ and

${\rm BR}(h_{2}\rightarrow XX)=\frac{\Gamma(h_{2}\rightarrow
XX)}{\sin^{2}\alpha\,{\Gamma^{\rm
SM}_{h}}|_{m_{h}=m_{h_{2}}}+\Gamma_{h_{2}}^{\rm DS}}\;.$ (48)

The partial decay with of $h_{2}$ to the dark sector particles
$\Gamma_{h_{2}}^{\rm DS}\simeq\Gamma(h_{2}\rightarrow
Z_{2}Z_{2})+\Gamma(h_{2}\rightarrow Z_{3}Z_{3})+\Gamma(h_{2}\rightarrow{\cal
W}^{p}{\cal W}^{m})$. We utilize the current direct heavy resonance searches
at the LHC including $h_{2}\rightarrow W^{+}W^{-},ZZ$ ATLAS:2020fry ,
$h_{2}\rightarrow\tau^{+}\tau^{-}$ ATLAS:2020zms  and
$h_{2}\rightarrow\mu^{+}\mu^{-}$ ATLAS:2019odt  to set constraints on the
model parameter space. We find that the current data from $\tau^{+}\tau^{-}$
and $\mu^{+}\mu^{-}$ final states are unlikely to constraint the model,
however the di-boson final states data can put an upper bound on the mixing
angle as depicted as purple and green shaded regions in Fig. 1. We assume
$\Gamma_{h_{2}}^{\rm DS}\sim 0$ for the result shown in Fig. 1, and notice
that a larger value of this partial decay width leads to a less stringent
upper bound on the mixing angle from the heavy Higgs boson searches.

###  III.2 Electroweak Precision Measurements and Dark Photon/$Z^{\prime}$
Physics

The presence of the extra scalar $h_{2}$ and its mixing to the SM-like Higgs
boson can modify the oblique parameters $S$, $T$, and $U$ Peskin:1991sw . In
particular, any oblique parameter $\cal O$ in the model can be given as
Profumo:2014opa

$\Delta{\cal O}=\left[{\cal O}^{\rm SM}(m_{h_{2}})-{\cal O}^{\rm
SM}(m_{h_{1}})\right]\,\sin^{2}\alpha\;,$ (49)

where ${\cal O}^{\rm SM}$ is the oblique parameter given in the SM. We utilize
the global fit values for the oblique parameters taken from the Particle Data
Group (PDG) Zyla:2020zbs , which are given as

$\displaystyle\Delta S$ $\displaystyle=$ $\displaystyle-0.01\pm 0.1\;,$
$\displaystyle\Delta T$ $\displaystyle=$ $\displaystyle 0.03\pm 0.12\;,$ (50)
$\displaystyle\Delta U$ $\displaystyle=$ $\displaystyle 0.02\pm 0.11\;,$

and the correlation coefficients are $0.92,-0.8$ and $-0.93$ for ($\Delta
S,\Delta T$), ($\Delta S,\Delta U$) and ($\Delta T,\Delta U$), respectively.
The constraint on the mixing angle $\alpha$ versus the heavy Higgs mass from
the oblique parameters are shown as the blue shaded region in Fig. 1. One can
see that the upper bound on the mixing angle becomes more stringent in the
heavier mass region of $h_{2}$. At $m_{h_{2}}=1$ TeV, the mixing angle
$|\sin\alpha|\lesssim 0.23$.

The presence of extra neutral gauge bosons and its mixing with the SM gauge
boson given in (39) can modify the physics at $Z$-pole. The most relevant
parameter to the modification is the kinetic mixing $\epsilon$. We compute the
$Z$-pole physics shifts and use the electroweak precision measurements given
in PDG Zyla:2020zbs  to constraint the model.

Moreover, via the kinetic mixing, the dark photon and dark $Z^{\prime}$ can
decay to the SM particles if they are kinematically allowed. Their decay
widths to the SM leptons can be given as

$\Gamma\left(Z_{i}\rightarrow\bar{\ell}\ell\right)=\frac{\alpha}{3}\varepsilon_{\ell}^{2}m_{Z_{i}}\sqrt{1-\mu_{\ell}^{2}}\left(1+\frac{\mu_{\ell}^{2}}{2}\right),$
(51)

where $i=\\{2,3\\}$, $\mu_{\ell}=2m_{\ell}/m_{Z_{i}}<1$ since it only opens
for $m_{Z_{i}}>2m_{\ell}$ and the mixing parameter $\varepsilon_{\ell}$ is
given at tree level by

$\varepsilon_{\ell}=\frac{1}{2gs_{W}}\sqrt{\left(C_{V}^{\ell
i}\right)^{2}+\left(C_{A}^{\ell
i}\right)^{2}\left(\frac{1-\mu_{\ell}^{2}}{1+\mu_{\ell}^{2}/2}\right)}\;,$
(52)

with $C_{V}^{\ell i}$ and $C_{A}^{\ell i}$ are given in (87) and (88),
respectively. Current searches for the dark photon and dark $Z^{\prime}$ set
constraints on the mixing parameter $\varepsilon_{\ell}$ and $m_{Z_{2,3}}$. We
focus on the mass regions $m_{Z_{2,3}}>1$ MeV and take into account the
current dark photon searches from collider/fixed target experiments including
A1 Merkel:2014avp , KLOE KLOE-2:2011hhj ; KLOE-2:2012lii ; KLOE-2:2014qxg ;
KLOE-2:2016ydq , NA48/2 NA482:2015wmo , BaBar BaBar:2014zli , LHCb
LHCb:2019vmc , CMS CMS:2019buh  and ATLAS ATLAS:2019erb , and beam dump
experiments including E774 Bross:1989mp , E141 Riordan:1987aw , E137
Bjorken:1988as ; Batell:2014mga ; Marsicano:2018krp ), $\nu$-Cal
Blumlein:2011mv ; Blumlein:2013cua , and CHARM Gninenko:2012eq , as well as
the bounds from supernovae Chang:2016ntp  and $(g-2)_{e}$ Pospelov:2008zw .
Recent constraints from FASER FASER:2023tle  (decaying to $e^{+}e^{-}$) and
NA62 NA62:2023qyn  (decaying to $\mu^{+}\mu^{-}$) experiments are also taken
into account in this analysis.

###  III.3 Relic Density

The DM candidate ${\cal W}$ can annihilate to SM particles via the
$2\rightarrow 2$ $s$ channel processes with $h_{i}$ and $Z_{i}$ as mediators.
In addition, they can annihilate to $Z_{i}Z_{j}$ and $h_{i}h_{j}$ via seagull,
$t$ channel and $u$ channel diagrams.

Here we utilize micrOMEGAs package Belanger:2018ccd  to calculate the DM relic
density and compare against the observed DM relic density, $\Omega_{\rm
DM}h^{2}=0.120\pm 0.001$ from the Planck Collaboration Planck:2018vyg .

We note that the annihilation process ${\cal W}^{p}{\cal
W}^{m}\rightarrow\gamma\gamma$ can be occurred via the SM charged fermion and
$W^{\pm}$ loops (triangle diagrams with $h_{i}/Z_{j}$ exchange connected to
the $\mathcal{W}$s) and therefore negligible in the relic abundance
computation! It may give rise to discrete gamma ray lines in the sky for
indirect detection which is important and deserves a separate study by itself.

###  III.4 Direct Detection

The DM candidate ${\cal W}$ can be directly observable when it scatters off
the target material within underground detectors, resulting in a detectable
recoil energy. The scattering process between ${\cal W}$ and nucleons can
occur through $t$ channel diagrams involving mediators $h_{i}$ and $Z_{i}$ in
the model. The elastic cross section between DM and a nucleon, denoted as $N$,
can be expressed as

$\displaystyle\sigma^{\rm SI}_{{\cal W}N}$ $\displaystyle=\sigma^{\rm
SI}_{{\cal
W}p}\frac{\sum_{k}\eta_{k}\mu_{A_{k}}^{2}\left[Z_{\text{atom}}+(A_{k}-Z_{\text{atom}})f_{n}/f_{p}\right]^{2}}{\sum_{k}\eta_{k}\mu_{A_{k}}^{2}A_{k}^{2}}\;.$
(53)

Here $\sigma^{\rm SI}_{{\cal W}p}$ represents the elastic cross section
between DM and a proton, $\mu_{p}$ is the reduced DM-proton mass defined as
$\mu_{p}=m_{{\cal W}}m_{p}/(m_{{\cal W}}+m_{p})$, $\mu_{A_{k}}=m_{{\cal
W}}m_{A_{k}}/(m_{{\cal W}}+m_{A_{k}})$ is the reduced DM-isotope nucleus mass
and $f_{p}$ and $f_{n}$ are effective couplings of the DM with protons and
neutrons, respectively. The atomic number is denoted as $Z_{\text{atom}}$,
while $\eta_{k}$ and $A_{k}$ represent isotope-specific variables indicating
the abundance and mass number of the $k^{\text{th}}$ target isotope,
respectively.

Typically, direct detection experiments assume isospin conservation, i.e.,
$f_{p}=f_{n}$, and report results in the form of (53). However, in this model,
the couplings between the gauge boson mediator particles $Z_{i}$ and quarks,
specifically $u$ and $d$, can differ due to their distinct SM charges, leading
to isospin violation, where $f_{p}\neq f_{n}$. Following Refs. Feng:2011vu ;
Yaguna:2016bga , we can rescale the reported experimental limit,
$\sigma_{\text{limit}}\rightarrow\sigma_{\text{limit}}\times\sigma^{\rm
SI}_{{\cal W}p}/\sigma^{\rm SI}_{{\cal W}N}$ to account for the effects of
isospin violation. This rescaled value can be used to constrain $\sigma^{\rm
SI}_{{\cal W}p}$.

We use micrOMEGAs package to calculate the effective couplings $f_{p}$ and
$f_{n}$ as well as the elastic cross section $\sigma^{\rm SI}_{{\cal W}p}$.
The results of the cross section are compared against recent upper limits from
experiments including CRESST III CRESST:2017ues , DarkSide-50 DarkSide:2018bpj
; DarkSide-50:2022qzh , XENON1T XENON:2018voc ; XENON:2019gfn , XENONnT
XENON:2023cxc , PandaX-4T PandaX-4T:2021bab , and LZ LZ:2022lsv .

Figure 2:  The seagull diagram for SIDM in the model. Figure 3:  The $t$
channel Feynman diagrams for SIDM in the model.

##  IV Self Interacting Dark Matter

The non-Abelian nature of gauge bosons enables the self-interaction of DM
$\mathcal{W}$ via the seagull diagram as depicted in Fig. 2. Furthermore,
$\mathcal{W}$ self-interaction can be facilitated through $t$ channel diagrams
involving the $h_{i}$ and $Z_{i}$ bosons as mediators, depicted in Fig. 3.
While the seagull diagram interaction results in a velocity-independent cross
section, the mediator exchange diagrams can give rise to long-range
interactions with velocity-dependence. When the mediator involves a vector
boson (i.e., $Z_{i}$ exchange), interactions between the DM particle ${\cal
W}^{m}$ and its antiparticle ${\cal W}^{p}$ result in an attractive potential,
while interactions between two identical particles ${\mathcal{W}}^{m}$ (or two
identical antiparticles ${\cal W}^{p}$) yield a repulsive potential.
Conversely, if the mediator is a scalar boson (i.e., $h_{i}$ exchange), the DM
self-interactions manifest as purely attractive. Since we are considering a
symmetric DM (equal ${\cal W}^{m}$ and ${\cal W}^{p}$), it exhibits both
attractive and repulsive self-interactions DM and hence the total effective
cross section will be given as the average of the two.

The momentum transfer cross section for the DM self interacting is obtained by

$\sigma_{\rm T}\equiv\int d\Omega\,(1-\cos\theta)\frac{d\sigma}{d\Omega}\;,$
(54)

where $\theta$ is the scattering angle and $(1-\cos\theta)$ is the fractional
change in longitudinal momentum. The cross section is calculated according to
different regimes defined by two dimensionless parameters $\kappa$ and $\beta$
introduced in (4) and presented again here for convenience

$\kappa=\frac{m_{\cal
W}v}{m_{X}}\;,\quad\quad\beta=\frac{2\alpha_{X}m_{X}}{m_{\cal W}v^{2}}\;,$
(55)

where $X\equiv\\{{h_{i}},{Z_{i}}\\}$ in the model and thus

$\alpha_{h_{1}}=\frac{1}{4\pi}\left(\frac{g_{D}\sin\alpha}{2}\right)^{2}\,,\;\;\;\;\alpha_{h_{2}}=\frac{1}{4\pi}\left(\frac{g_{D}\cos\alpha}{2}\right)^{2}\,,\;\;\;\;\alpha_{Z_{i}}=\frac{\left(g_{D}{\cal
O}_{3i}\right)^{2}}{4\pi}\,.$ (56)

In particular, the momentum transfer cross section can be computed in the Born
approximation regime if $2\beta\kappa^{2}\ll 1$, in the semi-classical regime
if $\kappa\geq 1$ and in the quantum regime if $\kappa\ll 1$.

In the Born approximation, the cross section can be computed perturbatively in
$\alpha_{X}$. On the other hand, in the semi-classical and quantum regimes,
the non-perturbative effects become significant. We follow the calculation
method in Colquhoun:2020adl  for the semi-classical regime. In particular, the
momentum transfer cross section in the semi-classical regime for an attractive
Yukawa potential is given by

$\displaystyle\sigma_{\mathrm{T}}^{\mathrm{att.}}$
$\displaystyle=\frac{\pi}{m_{X}^{2}}\times\begin{cases}2\beta^{2}\zeta_{1/2}\left(\kappa,\beta\right)&\beta\leq
0.2\,,\\\ \hskip 199.16928pt\ &\\\\[-11.38109pt]
2\beta^{2}\zeta_{1/2}\left(\kappa,\beta\right)e^{0.64(\beta-0.2)}&0.2<\beta\leq
1\,,\\\ \hskip 199.16928pt\ &\\\\[-11.38109pt]
4.7\log(\beta+0.82)&1<\beta<50\,,\\\ \hskip 199.16928pt\ &\\\\[-11.38109pt]
2(\log\beta)(\log\log\beta+1)&\beta\geq 50\;,\end{cases}$ (57)

and for a repulsive Yukawa potential,

$\displaystyle\sigma_{\mathrm{T}}^{\mathrm{rep.}}$
$\displaystyle=\frac{\pi}{m_{X}^{2}}\times\begin{cases}2\beta^{2}\zeta_{1/2}\left(\kappa,\beta\right)&\beta\leq
0.2\,,\\\ \hskip 199.16928pt\ &\\\\[-11.38109pt]
2\beta^{2}\zeta_{1/2}\left(\kappa,\beta\right)e^{-0.53(\beta-0.2)}&0.2<\beta\leq
1\,,\\\ \hskip 199.16928pt\ &\\\\[-11.38109pt]
2.9\log(\beta+0.47)&1<\beta<50\,,\\\ \hskip 199.16928pt\ &\\\\[-11.38109pt]
\lambda_{\mathrm{T}}(\log 2\beta-\log\log 2\beta)^{2}&\beta\geq
50\;,\end{cases}$ (58)

where

$\displaystyle\lambda_{\mathrm{T}}$ $\displaystyle=\frac{1}{2}(1+\cos 2+2\sin
2)\;,$ (59) $\displaystyle\zeta_{n}(\kappa,\beta)$
$\displaystyle=\frac{\text{max}(n,\beta\kappa)^{2}-n^{2}}{2\kappa^{2}\beta^{2}}+\eta\left(\frac{\text{max}(n,\beta\kappa)}{\kappa}\right)\;,$
(60) $\displaystyle\eta(x)$
$\displaystyle=x^{2}\left[-K_{1}\left(x\right)^{2}+K_{0}\left(x\right)K_{2}\left(x\right)\right]\;,$
(61)

with $n=1$ for distinguishable particles and $n=\tfrac{1}{2}$
($n=\tfrac{3}{2}$) for identical particles with even (odd) spatial wave
function, and $K_{n}(x)$ being the modified Bessel function of the second kind
of order $n$.

For the quantum regime where $\kappa<0.4$, the cross section is largely
dominated by $S$-wave scattering and the Schrödinger equation can be solved
analytically for the Hulthén potential Tulin:2013teo ,

$\sigma_{T}^{\mbox{Hulth\'{e}n}}=\frac{16\pi}{m_{\cal
W}^{2}v^{2}}\sin^{2}\delta_{0}\;,$ (62)

where the phase shift $\delta_{0}$ is given by

$\delta_{0}=\arg\left(i\frac{\Gamma(l_{p}+l_{m}-2)}{\Gamma(l_{p})\Gamma(l_{m})}\right)\;,$
(63)

with

$\begin{split}l_{p}&=1+i\frac{\kappa}{{2\rho}}(1+\sqrt{1\pm 2\rho\beta})\;,\\\
l_{m}&=1+i\frac{\kappa}{{2\rho}}(1-\sqrt{1\pm 2\rho\beta})\;.\end{split}$ (64)

Here $\rho\simeq 1.6$ is a dimensionless number and the signs $+(-)$ denotes
repulsive (attractive).

For the region between semi-classical and quantum regimes, i.e.
$0.4<\kappa<1$, we follow Colquhoun:2020adl  to use the interpolation function

$\displaystyle\sigma_{T}=(1-\kappa)\,\sigma_{T}^{\mbox{Hulth\'{e}n}}/0.6+(\kappa-0.4)\,\sigma_{T}^{\rm{rep(att)}}/0.6\;.$
(65)

Finally, to compare against the astrophysical data, we compute the velocity-
averaged transfer cross section as Colquhoun:2020adl

$\overline{\sigma_{T}}=\frac{\langle\sigma_{T}v^{2}\rangle}{16\sqrt{2}v_{0}^{2}/\pi}=\frac{1}{16\sqrt{2}v_{0}^{2}/\pi}\int_{0}^{\infty}f(v,\sqrt{2}v_{0})\sigma_{T}v^{2}dv\;,$
(66)

where $v$ is the relative velocity in the center-of-mass frame,
$v_{0}=\sqrt{\pi}\langle v\rangle/4$ is the velocity dispersion and the
Maxwell-Boltzmann distribution with parameter $a$ is given by

$f(v,a)=\sqrt{\frac{2}{\pi}}\frac{v^{2}e^{-v^{2}/(2a^{2})}}{a^{3}}\;.$ (67)

Figure 4: The momentum transfer cross section divided the DM mass from the
seagull diagram spanned in ($m_{\cal W},g_{D}$) plane. The blue-shaded region
is favored for the astrophysical observations in dwarfs, LSB spiral galaxies
and galaxy clusters. The dashed green, blue and red lines represent
$\sigma_{T}/{m_{\cal W}}=0.1,1$ and $10\,\rm{cm}^{2}/g$ respectively.

##  V Numerical Results

We begin by showing in Fig. 4 the Born cross section solely from the seagull
diagram spanned in the ($m_{\cal W},g_{D}$) plane. The blue band indicates the
parameter space where $\sigma_{T}/m_{\cal W}=0.1-10\,\rm{cm}^{2}/g$, as
preferred for solving small scale issues. The dark matter mass $m_{\cal W}$
lies in sub-GeV scale for the new gauge coupling $g_{D}$ of order unity.

Next, we explore the parameter space of the model that satisfies all the
theoretical and experimental constraints discussed in the previous sections.
In particular, we include the theoretical constraints on the scalar potential,
the Higgs measurements, the direct searches for heavy Higgs, the electroweak
precision measurements, dark photon, dark $Z^{\prime}$ searches and the dark
matter relic density as well as the dark matter direct detection. We compute
the self interacting DM cross section as discussed in section IV and then
focus on the parameter space that give
$0.1\,{\rm{cm}^{2}/g}<\sigma_{T}/m_{\cal W}\,(\langle v\rangle=10\,{\rm
km/s})<10\,\rm{cm}^{2}/g$ as preferred for solving the small scale issues. To
sample the parameter space in the model, we employ MCMC scans using emcee
Foreman-Mackey:2012any . The parameter range is set as follows,

$\displaystyle m_{h_{2}}/{\rm GeV}$ $\displaystyle\in$
$\displaystyle[150,\,2000]\;,$ (68) $\displaystyle\alpha/{\rm rad}$
$\displaystyle\in$ $\displaystyle\left[-0.2,\,0.2\right]\;,$ (69)
$\displaystyle m_{{\cal W}}/{\rm GeV}$ $\displaystyle\in$
$\displaystyle[10^{-3},\,10^{3}]\;,$ (70) $\displaystyle M_{X}/{\rm GeV}$
$\displaystyle\in$ $\displaystyle[10^{-3},\,80]\;,$ (71) $\displaystyle g_{D}$
$\displaystyle\in$ $\displaystyle[10^{-6},\,1]\;,$ (72) $\displaystyle
g_{D}^{\prime}$ $\displaystyle\in$ $\displaystyle[10^{-6},\,1]\;,$ (73)
$\displaystyle\epsilon$ $\displaystyle\in$
$\displaystyle[10^{-15},\,10^{-1}]\;.$ (74)

Here $m_{h_{2}}$ and $\alpha$ are sampled in the linear scale, while the
remaining free parameters are in the log scale prior. With the above choice of
scan parameters, $v_{D}$ is given by $2m_{\mathcal{W}}/g_{D}$, while
$\lambda$, $\lambda_{D}$ and $\lambda^{\prime}$ are computed by using (26),
(27) and (28) respectively. We note that the lower scanning range for
$\epsilon$ is set to a very small value due to the stringent constraint
imposed on the dark photon within the low mass region of $Z_{2,3}$, where
$m_{Z_{2,3}}\sim{\cal O}(100\,{\rm MeV})$, as per the E137 Bjorken:1988as ;
Batell:2014mga ; Marsicano:2018krp , $\nu$-Cal Blumlein:2011mv ;
Blumlein:2013cua  and CHARM Gninenko:2012eq  experiments as well as the
supernova bound Chang:2016ntp .

Figure 5: The $2\sigma$ favored region projected on the ($m_{Z_{2}}$,
$m_{Z_{3}}$) plane. The blue boxes indicate the favored region considering all
model constraints except the SIDM, while the red and green stars represent the
inclusion of SIDM with the dominance of $Z_{2}$ and $Z_{3}$ mediation
diagrams, respectively.

Fig. 5 shows the $2\sigma$ favored region ($\Delta\chi^{2}<5.99$) projected on
($m_{Z_{2}}$, $m_{Z_{3}}$) plane. To understand the impact of the SIDM, we
illustrate two superimposed regions: one taking into account the SIDM
constraint (regions marked by red and green stars) and another that does not
include the SIDM constraint (blue boxed region). We find that the dominant
contribution to the SIDM cross section arises from either the $Z_{2}$ or
$Z_{3}$ mediation diagram, depending on their respective mass spectra, with
the lighter mass one significantly influencing the magnitude of this
contribution. Upon scrutinizing the constraints imposed by astrophysical data
on SIDM cross section, discrete and isolated regions are discernible within
the parameter space, each indicative of distinctive primary contributions. In
particular, the region marked by red stars, where $m_{Z_{2}}$ < 400 MeV and
$m_{Z_{3}}>1$ GeV, corresponds to the $Z_{2}$ mediation diagram dominant
region while the green stars region where $m_{Z_{2}}>1$ GeV and
$m_{Z_{3}}<600$ MeV, represents the region that the $Z_{3}$ mediation diagram
dominant. Following we will show the results for the mass order relation
$m_{Z_{2}}<m_{Z_{3}}$, thereby focusing exclusively on the region where the
cross section is dominated by the $Z_{2}$ mediation diagram. The results for
the reverse mass order relation $m_{Z_{2}}>m_{Z_{3}}$ can be found in the
appendix B.

Figure 6: The $2\sigma$ favored region projected on the ($m_{\cal W}$,
$g_{D}$) plane (left panel) and the ($m_{Z_{2}}$, $g_{D}^{\prime}$) plane
(right panel). The blue boxes indicate the favored regions considering all
model constraints except the SIDM, while the red stars represent the inclusion
of SIDM.

Fig. 6 shows the favored region projected on 2D parameter space of the model.
Here, we again illustrate two superimposed regions as in Fig. 5 but by
assuming $m_{Z_{2}}<m_{Z_{3}}$, thus the regions marked by green stars
disappeared. In the left panel of Fig. 6, the favored region is projected on
the ($m_{\cal W}$, $g_{D}$) plane. The constraints observed on the upper and
lower bounds of the coupling parameter $g_{D}$ at a certain DM mass primarily
arise from the necessity to meet the DM relic density criteria. A larger value
of $g_{D}$ results in an increased DM annihilation cross section, consequently
leading to an underabundance of DM relic, while a decreased value (except for
the DM annihilation resonance regions) causes an overabundance of DM relics.
Notably, the discernible dip occurring around $m_{\cal W}\simeq 63$ GeV can be
attributed to the resonance facilitated by the SM-like Higgs in the DM
annihilation process. Additionally, resonance processes involving the heavy
Higgs $h_{2}$ are evident at $m_{\cal W}\geq 75$ GeV, particularly within the
realm of small $g_{D}$ couplings. We note that since we scan the $h_{2}$ mass
in the heavy region, the DM annihilation resonance process via $h_{2}$
exchange does not occur in lower DM mass regions.

Furthermore, it becomes apparent that the constraint imposed by SIDM
establishes not only upper and lower bounds for $g_{D}$ but also imposes a
lower bound on the DM mass itself. Specifically, in the absence of the SIDM
constraint, the viable range for the coupling $g_{D}$ extends in $[0.01,0.8]$
and $m_{\cal W}>0.4$ GeV. However, upon the inclusion of the SIDM constraint,
the range for $g_{D}$ falls within $[0.02,0.6]$ and $m_{\cal W}>0.7$ GeV.

In the right panel of Fig. 6, the favored region is projected on the
($m_{Z_{2}}$, $g_{D}^{\prime}$) plane. The influence of the SIDM constraint on
this parameter space is of considerable significance, demanding a notably
substantial $U(1)_{X}$ gauge coupling, $g_{D}^{\prime}$ and the mass region
for $Z_{2}$. Specifically, the imposition of the SIDM constraint leads to the
requirement of a relatively elevated range for $g_{D}^{\prime}$, falling
between $3\times 10^{-3}$ and $0.4$, and constrains the $m_{Z_{2}}$ values to
below $400$ MeV. This is due to the dominance of the SIDM cross section by the
$Z_{2}$ mediated process. To satisfy a large cross section preferred for
solving small scale issues, this process necessitates a relatively sizable
dark fine-structure coupling, $\alpha_{Z_{2}}$ (as defined in (56)), and
presumes a discernible mass hierarchy between the dark matter particle and its
mediator as previously studied in literature Tulin:2013teo ;
Kaplinghat:2015aga .

We note that $Z_{2}$ can play the role of a dark photon that couples to the SM
particles through kinetic mixing. Within the $Z_{2}$ mass allowing by the SIDM
constraint (i.e. $m_{Z_{2}}<400$ MeV), the constraints derived from beam-dump
experiments and supernova observations require the kinetic mixing parameter to
remain diminutive, e.g. $\epsilon\lesssim 10^{-7}$ at $m_{Z_{2}}=300$ MeV.
Consequently, this results in a notably weak interconnection between the dark
and SM sectors facilitated primarily through exchanges involving the neutral
gauge boson $Z_{i}$.

Figure 7: The SIDM cross section as a function of the average velocity. The
red, blue, and green data points represent the astrophysical data from dwarf
galaxies, LSB spiral galaxies, and galaxy clusters, respectively. The solid
purple line corresponds to the best fit from our model parameter space to the
dwarf galaxies and LSB spiral galaxies data. The dashed gray lines indicate
the velocity-independent cross sections.

Fig. 7 illustrates the fitting of the SIDM cross section from the model to the
astrophysical data sourced from dwarf galaxies, LSB spiral galaxies and galaxy
clusters. The astrophysical data are taken from Kaplinghat:2015aga . The
parameter set in the model, which satisfies all current constraints and
generates the best fit curve (depicted as the solid purple line in Fig. 7) for
dwarf galaxies and LSB spiral galaxies, is characterized by specific values:
$m_{\cal W}=6.1$ GeV, $M_{\mathcal{X}}=6.3$ MeV, $g_{D}=0.0446$,
$g_{D}^{\prime}=0.0423$, $\epsilon=2.8\times 10^{-14}$, $\alpha=-0.034$
radian, and $m_{h_{2}}=790$ GeV. However, extending this fitting to
incorporate high-velocity data from galaxy clusters poses a challenge for the
model, particularly in explaining the DM abundance from thermal freeze-out.
The SIDM cross section required to fit the high-velocity data demands a
smaller dark matter mass and relatively large coupling $g_{D}$. Consequently,
this leads to a lower DM relic density in comparison to the observed data from
Planck satellite.

Figure 8: The scattering cross section between DM and proton as a function of
the DM mass. The color scheme corresponds to the encoding used in Fig. 6. The
dashed green, black, cyan and lime lines represent the future sensitivities
from CDEX CDEX:2017kys , XENONnT XENON:2020kmp , DarkSide-20k
DarkSide-20k:2017zyg , and DARWIN DARWIN:2016hyl  detectors. The orange region
indicates the neutrino floor area.

Fig. 8 shows the favored region projected on the DM-proton spin-independent
cross section $\sigma^{\rm SI}_{{\cal W}p}$ and DM mass plane. We find that
the contributions to the DM-proton scattering cross section are mainly from
the diagrams mediated by the Higgs boson $h_{1}$, and $Z_{2}$ depending on
their masses and couplings to both the DM and the SM particles. It is
noteworthy that although the cross section $\sigma^{\rm SI}_{{\cal W}p}$
mediated by $Z_{2}$ can be suppressed due to the smallness of the kinetic
mixing parameter $\epsilon$, its effectiveness can be enhanced within the
lower mass spectrum of $Z_{2}$. Present constraints derived from the direct
detection of DM impose an upper limit on $\sigma^{\rm SI}_{{\cal W}p}$.
Interestingly, it is important to observe that this limit can appear less
stringent in comparison to the one reported by experimental papers, owing to
the consideration of isospin violation effects in our analysis. Of particular
interest is the potential to explore the favored region through forthcoming DM
direct detection experiments such as CDEX CDEX:2017kys , XENONnT XENON:2020kmp
, DarkSide-20k DarkSide-20k:2017zyg  and DARWIN DARWIN:2016hyl . However, it
is crucial to note that a portion of this favored region aligns with the
neutrino floor area, where distinguishing dark matter signals from neutrino
signals becomes notably challenging. Moreover, to distinguish the favored
regions with or without the SIDM constraint, one has to go deep inside the
neutrino floor.

##  VI Conclusions

In this study, we have conducted a thorough analysis of a novel physics model
that extends the SM through the incorporation of a dark version of the SM
gauge-Higgs sector. This model introduces a new gauge group, $SU(2)_{D}\times
U(1)_{D}$, along with the scalar $\Phi_{D}$, a doublet under $SU(2)_{D}$ and a
singlet under the SM $SU(2)_{L}\times U(1)_{Y}$. The interactions and mass
spectra of particles in the model were derived, revealing the two mass mixing
effects between the dark scalar and the SM-like Higgs as well as between the
three neutral gauge bosons. In addition, kinetic mixing between the two
factors of abelian $U(1)$s and a Stueckelberg mass of the hidden $U(1)_{D}$
are included. These setup give rise to three different portals connecting the
SM and dark sectors. Notably, the complex vector gauge boson in the dark
sector, as demonstrated in this work, can be a vectorial self-interacting dark
matter candidate.

Our investigation of the model incorporated both theoretical and current
experimental constraints. Theoretical constraints included adherence to
perturbative limits for couplings, perturbative unitarity, and the bounded-
from-below condition on the scalar potential. Collider phenomenological
constraints for the scalar sector were explored using current SM 125 GeV Higgs
data and direct searches for heavy Higgs at the LHC. The analysis found that
the most stringent upper bound on the mixing angle $|\sin\alpha|$ is set by
the 125 GeV Higgs signal strength data from ATLAS, limiting it to less than
0.2. The current dark photon/dark $Z^{\prime}$ searches from collider and
beam-dump experiments along with astrophysical observations have also taken
into account in this analysis. These constraints limit the kinetic mixing
parameter to be tiny, especially in the dark photon/dark $Z^{\prime}$ mass
below $\cal{O}$(100) MeV. The electroweak precision measurements on the
oblique parameter and $Z$-pole physics were also taken into account properly.

We analyzed the vectorial dark matter phenomenology in the model. Dark matter,
through Higgs and neutral gauge boson portals, exhibits annihilation to SM
particles, and, interestingly, to particles within the dark sector, such as
dark photon/dark $Z^{\prime}$ or dark Higgs, which subsequently decay to SM
particles. We found that the latter annihilation channel dominates when the
final state particles in the dark sector are light. Moreover, the resonance
annihilation process can be occurred and become dominant when the dark matter
is about a half of the mediator mass. We scanned over the parameter space and
computed the thermal relic density of the vector dark matter. Comparing
against the observed data from the Planck Collaboration, a linear relationship
between the dark matter mass and the $SU(2)_{D}$ gauge coupling is established
in non-resonance annihilation regions. A larger dark matter mass requires a
correspondingly larger gauge coupling, and at specific mass values, the gauge
coupling must fall within a certain range to align with observed data on dark
matter relic density. The vectorial dark matter mass ranges from sub-GeV to
TeV scale, with the $SU(2)_{D}$ gauge coupling varying between $\sim 10^{-2}$
and $0.7$.

The spin-independent cross section for the scattering process between dark
matter and nucleons in underground detectors was computed, accounting for
isospin violation effects. The primary contributions to the scattering cross
section arise from Higgs boson and dark photon exchange diagrams. To assess
the validity of our model, we incorporated the latest constraints from dark
matter direct detection experiments, including CRESST-III, DarkSide-50,
XENON1T, XENONnT, PandaX-4T, and LZ.

Our investigation extended to the self-interaction of vectorial dark matter
within the model, with the aim of addressing the small scale “crisis” observed
in astrophysical phenomena. Self-interaction in the model can occur through
the seagull contact interaction and exchange diagrams mediated by dark Higgs
and neutral gauge bosons. Notably, with an $SU(2)_{D}$ gauge coupling of order
unity and a dark matter mass around sub-GeV, the former yields a cross section
that aligns with resolving the small scale issues. However, this contact
interaction results in a velocity-independent cross section, posing challenges
in explaining observational data from astrophysical objects spread over vast
scales.

On the other hand, the mediator exchange interactions can give rise a
velocity-dependent cross section. Our findings indicate that with dark
photon/dark $Z^{\prime}$ boson masses below ${\cal O}$(100) MeV and new gauge
couplings ranging from ${\cal O}(10^{-2})$ to ${\cal O}(10^{-1})$, the cross
section falls within the expected range to address the small scale “crisis”.
Model parameters were fitted to astrophysical observations for dwarf galaxies,
low surface brightness spiral galaxies, and galaxy clusters. While the model
aligns well with low-velocity data from dwarf and low surface brightness
spiral galaxies, reconciling it with high-velocity data from galaxy clusters
poses a challenge, particularly when requiring the vectorial dark matter to
achieve the correct relic density from Planck data.

Interestingly, good portion of the favorable region of parameter space
consistent with all constraints, whether the vectorial dark matter has self-
interaction or not, is above the neutrino floor and therefore can be probed by
upcoming dark matter direct detection experiments such as CDEX, XENONnT,
DarkSide-20k, and DARWIN. This underscores the potential for experimental
validation and refinement of our model in light of evolving observational
capabilities.

## Acknowledgments

We would like to thank Professor Hai-Bo Yu for useful correspondence and Dr.
Sudhakantha Girmohanta for enlightening discussions. Our special thanks go to
Professor Tim Tait for useful suggestions. This work was supported in part by
the NSTC grant No. 111-2112-M-001-035 (TCY) and the National Natural Science
Foundation of China, grant No. 19Z103010239 (VQT). VQT would like to thank the
Medium and High Energy Group at the Institute of Physics, Academia Sinica,
Taiwan for their hospitality during the progress of this work.

##  Appendix A Neutral Current Interactions

In this appendix, we present the relevant neutral current interaction
Lagrangians for both the fermions and gauge bosons in the model. Our
conventions follow Peskin and Schroeder. 444M. E. Peskin and D. V. Schroeder,
“An Introduction to quantum field theory,” Addison-Wesley, 1995, ISBN
978-0-201-50397-5.

###  A.1 Fermions

We can express the covariant derivative acting onto any field multiplet $\psi$
in the model in terms of the physical fields $A$ and $Z_{i}(i=1,2,3)$. Recall
that in general we have 555 Here $T^{\pm}=T^{1}\pm iT^{2}$ and $T^{3}$ (and
similarly $\mathcal{T}^{\pm}=\mathcal{T}^{1}\pm i\mathcal{T}^{2}$ and
$\mathcal{T}^{3}$) where $T^{i}$ (and $\mathcal{T}^{i}$) $(i=1,2,3)$ are the
representation of the generators of $SU(2)$ depending on the field multiplet
$\psi$; $Y$ and $X$ are respectively the $U(1)_{Y}$ hypercharge and $U(1)_{D}$
dark hypercharge of $\psi$.

$\displaystyle D_{\mu}=\partial_{\mu}$
$\displaystyle-i\frac{g}{2}\left(W^{+}_{\mu}T^{+}+W^{-}_{\mu}T^{-}\right)-igW_{\mu}^{3}T^{3}-ig^{\prime}YB_{\mu}$
$\displaystyle-i\frac{g_{D}}{2}\left({\mathcal{W}}^{p}_{\mu}{\mathcal{T}}^{p}+{\mathcal{W}}^{m}_{\mu}{\mathcal{T}}^{m}\right)-ig_{D}{\mathcal{W}}^{3}_{\mu}\mathcal{T}^{3}-ig_{D}^{\prime}X{\mathcal{X}}_{\mu}\;.$
(75)

After the three consecutive transformations of $K$, $O_{W}$ and $\mathcal{O}$
obtained above, the 4 original neutral fields $B$, $W^{3}$, $\mathcal{X}$ and
$\mathcal{W}^{3}$ transform into physical fields $A$, $Z_{1}$, $Z_{2}$ and
$Z_{3}$. Thus one recasts the covariant derivative in (A.1) into

$\displaystyle D_{\mu}=\partial_{\mu}$ $\displaystyle-
ieQA_{\mu}-i\frac{g}{2}\left(W^{+}_{\mu}T^{+}+W^{-}_{\mu}T^{-}\right)-i\frac{g_{D}}{2}\left({\mathcal{W}}^{p}_{\mu}{\mathcal{T}}^{p}+{\mathcal{W}}^{m}_{\mu}{\mathcal{T}}^{m}\right)$
$\displaystyle-i\sum_{i=1}^{3}\left[\frac{g}{c_{W}}\left(T^{3}-s_{W}^{2}Q\right)\mathcal{O}_{1i}-\left(g^{\prime}S_{\epsilon}Y-g^{\prime}_{D}C_{\epsilon}X\right)\mathcal{O}_{2i}+g_{D}\mathcal{T}^{3}\mathcal{O}_{3i}\right]Z_{i}\;,$
(76)

where $Q=T^{3}+Y$ is the familiar electric charge operator in SM in unit of
$e=gs_{W}=g^{\prime}c_{W}$. Using (A.1), one can work out the couplings of the
SM fermions $\psi_{L}$ and $\psi_{R}$ with all the gauge fields. Note that for
the SM fermions their $X$ and $\mathcal{T}^{i}(i=1,2,3)$ are absent!

Explicitly, the neutral current couplings for the SM fermions are given by

$\mathcal{L}_{\rm
NC}\supset\sum_{i}\sum_{f}\bar{f}\gamma_{\mu}\left(C^{fi}_{L}P_{L}+C^{fi}_{R}P_{R}\right)fZ^{\mu}_{i}\;,$
(77)

where the summation is over all SM fermions $f$ and the three neutral gauge
bosons $Z_{i}$. The couplings $C^{fi}_{L}$ and $C^{fi}_{R}$ are given by

$\displaystyle C^{ui}_{L}$
$\displaystyle=\frac{g}{c_{W}}\left(\frac{1}{2}-\frac{2}{3}s_{W}^{2}\right)\mathcal{O}_{1i}-\frac{1}{6}g^{\prime}S_{\delta}\mathcal{O}_{2i}\;,$
(78) $\displaystyle C^{ui}_{R}$
$\displaystyle=\frac{g}{c_{W}}\left(-\frac{2}{3}s_{W}^{2}\right)\mathcal{O}_{1i}-\frac{2}{3}g^{\prime}S_{\delta}\mathcal{O}_{2i}\;,$
(79) $\displaystyle C^{di}_{L}$
$\displaystyle=\frac{g}{c_{W}}\left(-\frac{1}{2}+\frac{1}{3}s_{W}^{2}\right)\mathcal{O}_{1i}-\frac{1}{6}g^{\prime}S_{\delta}\mathcal{O}_{2i}\;,$
(80) $\displaystyle C^{di}_{R}$
$\displaystyle=\frac{g}{c_{W}}\left(\frac{1}{3}s_{W}^{2}\right)\mathcal{O}_{1i}+\frac{1}{3}g^{\prime}S_{\delta}\mathcal{O}_{2i}\;,$
(81) $\displaystyle C^{\nu i}_{L}$
$\displaystyle=\frac{g}{c_{W}}\left(\frac{1}{2}\right)\mathcal{O}_{1i}+\frac{1}{2}g^{\prime}S_{\delta}\mathcal{O}_{2i}\;,$
(82) $\displaystyle C^{\nu i}_{R}$ $\displaystyle=\;0\;,$ (83) $\displaystyle
C^{ei}_{L}$
$\displaystyle=\frac{g}{c_{W}}\left(-\frac{1}{2}+s_{W}^{2}\right)\mathcal{O}_{1i}+\frac{1}{2}g^{\prime}S_{\delta}\mathcal{O}_{2i}\;,$
(84) $\displaystyle C^{ei}_{R}$
$\displaystyle=\frac{g}{c_{W}}s_{W}^{2}\mathcal{O}_{1i}+g^{\prime}S_{\delta}\mathcal{O}_{2i}\;.$
(85)

Often $\mathcal{L}_{\rm NC}$ in (77) is written in an alternative form in the
literature as

$\mathcal{L}_{\rm
NC}\supset\frac{1}{2}\sum_{i}\sum_{f}\bar{f}\gamma_{\mu}\left(C^{fi}_{V}-C^{fi}_{A}\gamma_{5}\right)fZ^{\mu}_{i}\;,$
(86)

with

$\displaystyle C^{fi}_{V}$ $\displaystyle=C^{fi}_{L}+C^{fi}_{R}\;,$ (87)
$\displaystyle C^{fi}_{A}$ $\displaystyle=C^{fi}_{L}-C^{fi}_{R}\;.$ (88)

The couplings presented here are used to obtain the isospin violation effects
in direct detection.

###  A.2 Gauge Bosons

The cubic coupling between $W^{3}$ and $W^{+}W^{-}$ in the SM is given by

$\displaystyle\mathcal{L}_{\rm NC}\supset ig$
$\displaystyle\biggl{[}\left(\partial_{\mu}W^{+}_{\nu}-\partial_{\nu}W^{+}_{\mu}\right)W^{-\mu}W^{3\,\nu}-\left(\partial_{\mu}W^{-}_{\nu}-\partial_{\nu}W^{-}_{\mu}\right)W^{+\mu}W^{3\,\nu}\biggr{.}$
$\displaystyle\biggl{.}+\frac{1}{2}\left(\partial_{\mu}W^{3}_{\nu}-\partial_{\nu}W^{3}_{\mu}\right)\left(W^{+\mu}W^{-\nu}-W^{-\mu}W^{+\nu}\right)\biggr{]}\;.$
(89)

Using
$W^{3}=s_{W}A+c_{W}W^{3\,\prime\prime}=s_{W}A+c_{W}\mathcal{O}_{1i}Z_{i}$, we
obtain the same SM $\gamma W^{+}W^{-}$ coupling

$\displaystyle\mathcal{L}_{\rm NC}\supset ie$
$\displaystyle\biggl{[}\left(\partial_{\mu}W^{+}_{\nu}-\partial_{\nu}W^{+}_{\mu}\right)W^{-\mu}A^{\nu}-\left(\partial_{\mu}W^{-}_{\nu}-\partial_{\nu}W^{-}_{\mu}\right)W^{+\mu}A^{\nu}\biggr{.}$
$\displaystyle\biggl{.}+\frac{1}{2}\left(\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}\right)\left(W^{+\mu}W^{-\nu}-W^{-\mu}W^{+\nu}\right)\biggr{]}\;,$
(90)

where we have used the well-known relation $e=gs_{W}$. The $Z_{i}W^{+}W^{-}$
is given by

$\displaystyle\mathcal{L}_{\rm NC}\supset ig\sum_{i=1}^{3}\mathcal{O}_{1i}$
$\displaystyle\biggl{[}\left(\partial_{\mu}W^{+}_{\nu}-\partial_{\nu}W^{+}_{\mu}\right)W^{-\mu}Z^{\nu}_{i}-\left(\partial_{\mu}W^{-}_{\nu}-\partial_{\nu}W^{-}_{\mu}\right)W^{+\mu}Z^{\nu}_{i}\biggr{.}$
$\displaystyle\biggl{.}+\frac{1}{2}\left(\partial_{\mu}Z_{i\,\nu}-\partial_{\nu}Z_{i\,\mu}\right)\left(W^{+\mu}W^{-\nu}-W^{-\mu}W^{+\nu}\right)\biggr{]}\;.$
(91)

Similarly, the cubic coupling between $Z_{i}$ and
$\mathcal{W}^{p}\mathcal{W}^{m}$ is given by

$\displaystyle\mathcal{L}_{\rm NC}\supset
ig_{D}\sum_{i=1}^{3}\mathcal{O}_{3i}$
$\displaystyle\biggl{[}\left(\partial_{\mu}\mathcal{W}^{p}_{\nu}-\partial_{\nu}\mathcal{W}^{p}_{\mu}\right)\mathcal{W}^{m\mu}Z^{\nu}_{i}-\left(\partial_{\mu}\mathcal{W}^{m}_{\nu}-\partial_{\nu}\mathcal{W}^{m}_{\mu}\right)\mathcal{W}^{p\mu}Z^{\nu}_{i}\biggr{.}$
$\displaystyle\biggl{.}+\frac{1}{2}\left(\partial_{\mu}Z_{i\nu}-\partial_{\nu}Z_{i\mu}\right)\left(\mathcal{W}^{p\mu}\mathcal{W}^{m\nu}-\mathcal{W}^{m\mu}\mathcal{W}^{p\nu}\right)\biggr{]}\;.$
(92)

Recall that the SM Feynman rule for the
$\gamma(p_{3},\lambda)W^{+}(p_{2},\nu)W^{-}(p_{1},\mu)$ vertex is

$-ie\left[\left(p_{1}-p_{2}\right)_{\lambda}g_{\mu\nu}+\left(p_{2}-p_{3}\right)_{\mu}g_{\nu\lambda}+\left(p_{3}-p_{1}\right)_{\nu}g_{\lambda\mu}\right]\;,$
(93)

with all momenta flowing into the vertex. Thus the Feynman rules for the
$Z_{i}W^{+}W^{-}$ and $Z_{i}\mathcal{W}^{p}\mathcal{W}^{m}$ vertices can be
deduced simply by substituting $e$ in (93) with $g\mathcal{O}_{1i}$ and
$g_{D}\mathcal{O}_{3i}$ respectively. The mixing matrix element
$\mathcal{O}_{3i}$ connects the DM $\mathcal{W}^{p,m}$ with the SM particles
via the $Z_{i}$ exchange. Note that there is no
$\gamma\mathcal{W}^{p}\mathcal{W}^{m}$ since $W^{3}$ doesn’t mix with
$\mathcal{W}^{3}$. $\mathcal{W}^{p,m}$ carry only dark charge but no electric
charge.

##  Appendix B The Reverse Mass Ordering Case of $m_{Z_{2}}>m_{Z_{3}}$

Figure 9: The $2\sigma$ favored region projected on ($m_{\cal W}$, $g_{D}$)
plane (left panel) and ($m_{Z_{3}}$, $g_{D}^{\prime}$) plane (right panel) for
the case $m_{Z_{2}}>m_{Z_{3}}$. The blue boxes indicate the favored region
considering all model constraints except the self-interacting DM, while the
green stars represent the inclusion of self-interacting DM.

Here we present the numerical result for the case where $m_{Z_{2}}>m_{Z_{3}}$.
For this case the SIDM cross section dominated by the $Z_{3}$ mediation
diagram. The favored region spanned on the ($m_{\cal W}$, $g_{D}$) and
($m_{Z_{3}}$, $g_{D}^{\prime}$) planes are shown in the left and right panels
of Fig. 9, respectively. In comparison with the results for the
$m_{Z_{2}}<m_{Z_{3}}$ scenario depicted in Fig. 6, it requires a stronger
upper bound on $g_{D}$ and a stronger lower bound on $g_{D}^{\prime}$ within
the SIDM favored region.

## References

  * (1) D. N. Spergel and P. J. Steinhardt, “Observational evidence for selfinteracting cold dark matter,” Phys. Rev. Lett. 84, 3760-3763 (2000) [arXiv:astro-ph/9909386 [astro-ph]].
  * (2) V. Silveira and A. Zee, “Scalar Phantoms,” Phys. Lett. B 161, 136-140 (1985).
  * (3) K. Cheung, Y. L. S. Tsai, P. Y. Tseng, T. C. Yuan and A. Zee, “Global Study of the Simplest Scalar Phantom Dark Matter Model,” JCAP 10, 042 (2012) [arXiv:1207.4930 [hep-ph]].
  * (4) D. E. Holz and A. Zee, “Collisional dark matter and scalar phantoms,” Phys. Lett. B 517, 239-242 (2001) [arXiv:hep-ph/0105284 [hep-ph]].
  * (5) S. Tulin, H. B. Yu and K. M. Zurek, “Beyond Collisionless Dark Matter: Particle Physics Dynamics for Dark Matter Halo Structure,” Phys. Rev. D 87, no.11, 115007 (2013) [arXiv:1302.3898 [hep-ph]].
  * (6) M. Kaplinghat, S. Tulin and H. B. Yu, “Dark Matter Halos as Particle Colliders: Unified Solution to Small-Scale Structure Puzzles from Dwarfs to Clusters,” Phys. Rev. Lett. 116, no.4, 041302 (2016) [arXiv:1508.03339 [astro-ph.CO]].
  * (7) S. Tulin and H. B. Yu, “Dark Matter Self-interactions and Small Scale Structure,” Phys. Rept. 730, 1-57 (2018) [arXiv:1705.02358 [hep-ph]].
  * (8) B. Colquhoun, S. Heeba, F. Kahlhoefer, L. Sagunski and S. Tulin, “Semiclassical regime for dark matter self-interactions,” Phys. Rev. D 103, no.3, 035006 (2021) [arXiv:2011.04679 [hep-ph]].
  * (9) A. Aboubrahim, W. Z. Feng, P. Nath and Z. Y. Wang, “Self-interacting hidden sector dark matter, small scale galaxy structure anomalies, and a dark force,” Phys. Rev. D 103, no.7, 075014 (2021) [arXiv:2008.00529 [hep-ph]].
  * (10) W. C. Huang, Y. L. S. Tsai and T. C. Yuan, “G2HDM : Gauged Two Higgs Doublet Model,” JHEP 1604, 019 (2016) [arXiv:1512.00229 [hep-ph]].
  * (11) R. Ramos, V. Tran and T. C. Yuan, “A Sub-GeV Low Mass Hidden Dark Sector of SU(2)H× U(1)X,” JHEP 11, 112 (2021) [arXiv:2109.03185 [hep-ph]].
  * (12) R. Ramos, Van Que Tran and T. C. Yuan, “Complementary searches of low mass non-Abelian vector dark matter, dark photon, and dark $Z^{\prime}$,” Phys. Rev. D 103, no.7, 075021 (2021) [arXiv:2101.07115 [hep-ph]].
  * (13) M. L. Graesser, I. M. Shoemaker and L. Vecchi, “Asymmetric WIMP dark matter,” JHEP 10, 110 (2011) [arXiv:1103.2771 [hep-ph]].
  * (14) H. Iminniyaz, M. Drees and X. Chen, “Relic Abundance of Asymmetric Dark Matter,” JCAP 07, 003 (2011) [arXiv:1104.5548 [hep-ph]].
  * (15) T. Lin, H. B. Yu and K. M. Zurek, “On Symmetric and Asymmetric Light Dark Matter,” Phys. Rev. D 85, 063503 (2012) [arXiv:1111.0293 [hep-ph]].
  * (16) K. M. Zurek, “Asymmetric Dark Matter: Theories, Signatures, and Constraints,” Phys. Rept. 537, 91-121 (2014) [arXiv:1308.0338 [hep-ph]].
  * (17) B. Kors and P. Nath, “A Stueckelberg extension of the standard model,” Phys. Lett. B 586, 366-372 (2004) [arXiv:hep-ph/0402047 [hep-ph]].
  * (18) B. Kors and P. Nath, “Aspects of the Stueckelberg extension,” JHEP 07, 069 (2005) [arXiv:hep-ph/0503208 [hep-ph]].
  * (19) K. Cheung and T. C. Yuan, “Hidden fermion as milli-charged dark matter in Stueckelberg Z- prime model,” JHEP 03, 120 (2007) [arXiv:hep-ph/0701107 [hep-ph]].
  * (20) B. Holdom, “Two U(1)’s and Epsilon Charge Shifts,” Phys. Lett. B 166, 196-198 (1986).
  * (21) D. Feldman, Z. Liu and P. Nath, “The Stueckelberg Z-prime Extension with Kinetic Mixing and Milli-Charged Dark Matter From the Hidden Sector,” Phys. Rev. D 75, 115001 (2007) [arXiv:hep-ph/0702123 [hep-ph]].
  * (22) A. Arhrib, R. Benbrik and N. Gaur, “$H\rightarrow\gamma\gamma$ in Inert Higgs Doublet Model,” Phys. Rev. D 85, 095021 (2012) [arXiv:1201.2644 [hep-ph]].
  * (23) A. Arhrib, Y. L. S. Tsai, Q. Yuan and T. C. Yuan, “An Updated Analysis of Inert Higgs Doublet Model in light of the Recent Results from LUX, PLANCK, AMS-02 and LHC,” JCAP 06, 030 (2014) [arXiv:1310.0358 [hep-ph]].
  * (24) G. C. Branco, P. M. Ferreira, L. Lavoura, M. N. Rebelo, M. Sher and J. P. Silva, “Theory and phenomenology of two-Higgs-doublet models,” Phys. Rept. 516, 1-102 (2012) [arXiv:1106.0034 [hep-ph]].
  * (25) R. L. Workman et al. (Particle Data Group), Prog. Theor. Exp. Phys. 2022, 083C01 (2022). For the online version, see https://pdg.lbl.gov/.
  * (26) C. T. Huang, R. Ramos, V. Tran, Y. L. S. Tsai and T. C. Yuan, “Consistency of Gauged Two Higgs Doublet Model: Gauge Sector,” JHEP 09, 048 (2019) [arXiv:1905.02396 [hep-ph]].
  * (27) G. Aad et al. [ATLAS], Nature 607 (2022) no.7917, 52-59 [erratum: Nature 612 (2022) no.7941, E24] [arXiv:2207.00092 [hep-ex]].
  * (28) G. Aad et al. [ATLAS], Eur. Phys. J. C 80 (2020) no.12, 1165 [arXiv:2004.14636 [hep-ex]].
  * (29) G. Aad et al. [ATLAS], Phys. Rev. Lett. 125 (2020) no.5, 051801 [arXiv:2002.12223 [hep-ex]].
  * (30) M. Aaboud et al. [ATLAS], JHEP 07 (2019), 117 [arXiv:1901.08144 [hep-ex]].
  * (31) M. E. Peskin and T. Takeuchi, “Estimation of oblique electroweak corrections,” Phys. Rev. D 46 (1992), 381-409.
  * (32) S. Profumo, M. J. Ramsey-Musolf, C. L. Wainwright and P. Winslow, “Singlet-catalyzed electroweak phase transitions and precision Higgs boson studies,” Phys. Rev. D 91 (2015) no.3, 035018 [arXiv:1407.5342 [hep-ph]].
  * (33) P. A. Zyla et al. [Particle Data Group], “Review of Particle Physics,” PTEP 2020, no.8, 083C01 (2020).
  * (34) H. Merkel, P. Achenbach, C. Ayerbe Gayoso, T. Beranek, J. Bericic, J. C. Bernauer, R. Böhm, D. Bosnar, L. Correa and L. Debenjak, et al. Phys. Rev. Lett. 112, no.22, 221802 (2014) [arXiv:1404.5502 [hep-ex]].
  * (35) F. Archilli et al. [KLOE-2], Phys. Lett. B 706, 251-255 (2012) [arXiv:1110.0411 [hep-ex]].
  * (36) D. Babusci et al. [KLOE-2], Phys. Lett. B 720, 111-115 (2013) [arXiv:1210.3927 [hep-ex]].
  * (37) D. Babusci et al. [KLOE-2], Phys. Lett. B 736, 459-464 (2014) [arXiv:1404.7772 [hep-ex]].
  * (38) A. Anastasi et al. [KLOE-2], Phys. Lett. B 757, 356-361 (2016) [arXiv:1603.06086 [hep-ex]].
  * (39) J. R. Batley et al. [NA48/2], Phys. Lett. B 746, 178-185 (2015) [arXiv:1504.00607 [hep-ex]].
  * (40) J. P. Lees et al. [BaBar], Phys. Rev. Lett. 113, no.20, 201801 (2014) [arXiv:1406.2980 [hep-ex]].
  * (41) R. Aaij et al. [LHCb], Phys. Rev. Lett. 124, no.4, 041801 (2020) [arXiv:1910.06926 [hep-ex]].
  * (42) A. M. Sirunyan et al. [CMS], Phys. Rev. Lett. 124, no.13, 131802 (2020) [arXiv:1912.04776 [hep-ex]].
  * (43) G. Aad et al. [ATLAS], Phys. Lett. B 796, 68-87 (2019) [arXiv:1903.06248 [hep-ex]].
  * (44) A. Bross, M. Crisler, S. H. Pordes, J. Volk, S. Errede and J. Wrbanek, Phys. Rev. Lett. 67, 2942-2945 (1991).
  * (45) E. M. Riordan, M. W. Krasny, K. Lang, P. De Barbaro, A. Bodek, S. Dasu, N. Varelas, X. Wang, R. G. Arnold and D. Benton, et al. Phys. Rev. Lett. 59, 755 (1987).
  * (46) J. D. Bjorken, S. Ecklund, W. R. Nelson, A. Abashian, C. Church, B. Lu, L. W. Mo, T. A. Nunamaker and P. Rassmann, “Search for Neutral Metastable Penetrating Particles Produced in the SLAC Beam Dump,” Phys. Rev. D 38, 3375 (1988).
  * (47) B. Batell, R. Essig and Z. Surujon, “Strong Constraints on Sub-GeV Dark Sectors from SLAC Beam Dump E137,” Phys. Rev. Lett. 113, no.17, 171802 (2014) [arXiv:1406.2698 [hep-ph]].
  * (48) L. Marsicano, M. Battaglieri, M. Bondi’, C. D. R. Carvajal, A. Celentano, M. De Napoli, R. De Vita, E. Nardi, M. Raggi and P. Valente, “Dark photon production through positron annihilation in beam-dump experiments,” Phys. Rev. D 98, no.1, 015031 (2018) [arXiv:1802.03794 [hep-ex]].
  * (49) J. Blumlein and J. Brunner, “New Exclusion Limits for Dark Gauge Forces from Beam-Dump Data,” Phys. Lett. B 701, 155-159 (2011) [arXiv:1104.2747 [hep-ex]].
  * (50) J. Blümlein and J. Brunner, “New Exclusion Limits on Dark Gauge Forces from Proton Bremsstrahlung in Beam-Dump Data,” Phys. Lett. B 731, 320-326 (2014) [arXiv:1311.3870 [hep-ph]].
  * (51) S. N. Gninenko, “Constraints on sub-GeV hidden sector gauge bosons from a search for heavy neutrino decays,” Phys. Lett. B 713, 244-248 (2012) [arXiv:1204.3583 [hep-ph]].
  * (52) J. H. Chang, R. Essig and S. D. McDermott, “Revisiting Supernova 1987A Constraints on Dark Photons,” JHEP 01, 107 (2017) [arXiv:1611.03864 [hep-ph]].
  * (53) M. Pospelov, “Secluded U(1) below the weak scale,” Phys. Rev. D 80, 095002 (2009) [arXiv:0811.1030 [hep-ph]].
  * (54) H. Abreu et al. [FASER], [arXiv:2308.05587 [hep-ex]].
  * (55) E. Cortina Gil et al. [NA62], JHEP 09 (2023), 035 [arXiv:2303.08666 [hep-ex]].
  * (56) G. Bélanger, F. Boudjema, A. Goudelis, A. Pukhov and B. Zaldivar, “micrOMEGAs5.0: Freeze-in,” Comput. Phys. Commun. 231, 173-186 (2018) [arXiv:1801.03509 [hep-ph]].
  * (57) N. Aghanim et al. [Planck], Astron. Astrophys. 641, A6 (2020) [erratum: Astron. Astrophys. 652, C4 (2021)] [arXiv:1807.06209 [astro-ph.CO]].
  * (58) J. L. Feng, J. Kumar, D. Marfatia and D. Sanford, “Isospin-Violating Dark Matter,” Phys. Lett. B 703, 124-127 (2011) [arXiv:1102.4331 [hep-ph]].
  * (59) C. E. Yaguna, “Isospin-violating dark matter in the light of recent data,” Phys. Rev. D 95, no.5, 055015 (2017) [arXiv:1610.08683 [hep-ph]].
  * (60) G. Angloher et al. [CRESST], Eur. Phys. J. C 77, no.9, 637 (2017) [arXiv:1707.06749 [astro-ph.CO]].
  * (61) P. Agnes et al. [DarkSide], Phys. Rev. Lett. 121, no.8, 081307 (2018) [arXiv:1802.06994 [astro-ph.HE]].
  * (62) P. Agnes et al. [DarkSide-50], Phys. Rev. D 107, no.6, 063001 (2023) [arXiv:2207.11966 [hep-ex]].
  * (63) E. Aprile et al. [XENON], Phys. Rev. Lett. 121, no.11, 111302 (2018) [arXiv:1805.12562 [astro-ph.CO]].
  * (64) E. Aprile et al. [XENON], Phys. Rev. Lett. 123, no.25, 251801 (2019) [arXiv:1907.11485 [hep-ex]].
  * (65) E. Aprile et al. [XENON], Phys. Rev. Lett. 131, no.4, 041003 (2023) [arXiv:2303.14729 [hep-ex]].
  * (66) Y. Meng et al. [PandaX-4T], Phys. Rev. Lett. 127, no.26, 261802 (2021) [arXiv:2107.13438 [hep-ex]].
  * (67) J. Aalbers et al. [LZ], Phys. Rev. Lett. 131, no.4, 041002 (2023) [arXiv:2207.03764 [hep-ex]].
  * (68) D. Foreman-Mackey, D. W. Hogg, D. Lang and J. Goodman, “emcee: The MCMC Hammer,” Publ. Astron. Soc. Pac. 125, 306-312 (2013) [arXiv:1202.3665 [astro-ph.IM]].
  * (69) H. Ma et al. [CDEX], J. Phys. Conf. Ser. 1342, no.1, 012067 (2020) [arXiv:1712.06046 [hep-ex]].
  * (70) E. Aprile et al. [XENON], JCAP 11, 031 (2020) [arXiv:2007.08796 [physics.ins-det]].
  * (71) C. E. Aalseth et al. [DarkSide-20k], Eur. Phys. J. Plus 133, 131 (2018) [arXiv:1707.08145 [physics.ins-det]].
  * (72) J. Aalbers et al. [DARWIN], JCAP 11, 017 (2016) [arXiv:1606.07001 [astro-ph.IM]].