# A Sub-GeV Low Mass Hidden Dark Sector of $SU(2)_{H}\times U(1)_{X}$

Raymundo Ramos<EMAIL_ADDRESS>Institute of Physics, Academia
Sinica, Nangang, Taipei 11529, Taiwan  Van Que Tran<EMAIL_ADDRESS>School
of Physics, Nanjing University, Nanjing 210093, China  Tzu-Chiang Yuan
<EMAIL_ADDRESS>Institute of Physics, Academia Sinica, Nangang,
Taipei 11529, Taiwan

###### Abstract

We present a detailed study of the non-abelian vector dark matter candidate
$W^{\prime}$ with a MeV–GeV low mass range, accompanied by a dark photon
$A^{\prime}$ and a dark $Z^{\prime}$ of similar masses, in the context of a
gauged two-Higgs-doublet model with the hidden gauge group that has the same
structure as the Standard Model electroweak gauge group. The stability of dark
matter is protected by an accidental discrete $Z_{2}$ symmetry ($h$-parity)
which was usually imposed ad hoc by hand. We examine the model by taking into
account various experimental constraints including dark photon searches at
NA48, NA64, E141, $\nu$-CAL, BaBar and LHCb experiments, electroweak precision
data from LEP, relic density from Planck satellite, direct (indirect)
detection of dark matter from CRESST-III, DarkSide-50, XENON1T (Fermi-LAT),
and collider physics from the LHC. The theoretical requirements of bounded
from below of the scalar potential and tree level perturbative unitarity of
the scalar sector are also imposed. The viable parameter space of the model
consistent with all the constraints is exhibited. While a dark $Z^{\prime}$
can be the dominant contribution in the relic density due to resonant
annihilation of dark matter, a dark photon is crucial to dark matter direct
detection. We also demonstrate that the parameter space can be further probed
by various sub-GeV direct dark matter experimental searches at CDEX, NEWS-G
and SuperCDMS in the near future.

##  I Introduction

The stability of heavy dark matter (DM) is usually implemented in many
particle physics models beyond the standard model (SM) by imposing a discrete
$Z_{2}$ symmetry in the Lagrangian. These models include the simplest scalar
phantom model by adding just a singlet $Z_{2}$-odd scalar field (real or
complex) to the SM Silveira:1985rk , the popular inert two-Higgs-doublet model
(I2HDM) Deshpande:1977rw ; Ma:2006km , the minimal supergravity standard model
with $R$-parity (MSSM) Chamseddine:1982jx ; Nath:1982zq ; Nilles:1983ge ,
little Higgs model with $T$-parity (LHM) Cheng:2003ju ; Cheng:2004yc ;
Low:2004xc  etc. In I2HDM, the DM candidate can be either the CP-even or -odd
scalar residing in the second $Z_{2}$-odd Higgs doublet. Many detailed
analysis of DM phenomenology in the scalar phantom models and I2HDM can be
found in the literature in McDonald:1993ex ; Burgess:2000yq ; Cheung:2012xb
and Barbieri:2006dq ; LopezHonorez:2006gr ; Arhrib:2013ela ; Belyaev:2016lok ;
Lu:2019lok ; Fabian:2020hny  respectively. In MSSM, the lightest
supersymmetric particle (LSP) for the DM candidate can be the spin 0 sneutrino
(the superpartner of neutrino) or the lightest spin 1/2 neutralino
Akula:2011aa  (in general a linear combination of two gauginos and two
Higgsinos Arnowitt:1995vg ). We note that in some low energy supergravity
models, the LSP can be the spin 3/2 gravitino, the superpartner of graviton.
For a review of supersymmetric dark matter, see for example Jungman:1995df .
In LHM, the spin 1 $T$-odd partner of the photon can be the DM candidate whose
collider implication was studied in Chen:2006ie . There are also well-
motivated non-abelian dark matter models based on additional gauge group like
$SU(2)$ Carone:2013wla ; Davoudiasl:2013jma ; Barman:2017yzr ; Barman:2019lvm
; Abe:2020mph ; Hu:2021pln , in which the extra spin 1 gauge boson
$W^{\prime}$ can be a DM candidate. Moreover, instead of specifying an
underlying dark matter model, one can also use the effective dark matter
theory approach Cao:2009uw ; Goodman:2010ku ; Cheung:2012gi  to discuss
various dark matter phenomenologies Goodman:2010yf ; Goodman:2010qn ;
Cheung:2010ua ; Cheung:2011nt ; Huang:2019ikw ; NathMaity:2021cne .

Recently a gauged two-Higgs-doublet model (G2HDM) based on an extended
electroweak gauge group $SU(2)_{L}\times U(1)_{Y}\times SU(2)_{H}\times
U(1)_{X}$, was proposed Huang:2015wts , in which a hidden discrete $Z_{2}$
symmetry ($h$-parity) Chen:2019pnt  arises naturally as an accidental symmetry
rather than imposed by hand. This discrete symmetry ensures the stability of
the DM candidate in G2HDM, which can be either a complex scalar (which in
general is a linear combination of various fields in the model) or a heavy
neutrino $\nu^{H}$ or an extra gauge boson $W^{\prime(p,m)}$, all of which
have odd $h$-parity. We note that, unlike the left-right symmetric model
Mohapatra:1979ia ; Mohapatra:1980yp , the $W^{\prime(p,m)}$ in G2HDM do not
carry electric charge.

The novel idea of G2HDM, as compared with many variants of general 2HDM
Branco:2011iw , is that the two Higgs doublets $H_{1}$ and $H_{2}$ of
$SU(2)_{L}$ are grouped into a fundamental representation of a new gauge group
$SU(2)_{H}$. Consistency checks of the model were scrutinized for the scalar
and gauge sectors in Arhrib:2018sbz  and Huang:2019obt  respectively. In
Chen:2019pnt , a detailed phenomenological analysis of the scalar DM candidate
in G2HDM was carried out. In general, the scalar DM candidate is a complex
field made up of a linear combination of the components from the inert Higgs
doublet $H_{2}$ and the $SU(2)_{H}$ doublet $\Phi_{H}$ and triplet
$\Delta_{H}$. By performing a detailed parameter scan in the model it was
demonstrated Chen:2019pnt  that only the triplet-like DM is favored when all
the constraints from the relic density, direct and indirect searches for the
DM are taken into account. As discussed in Huang:2015wts , the triplet
$\Delta_{H}$ plays the primary role as a trigger for spontaneous symmetry
breaking of the model down to $U(1)_{\rm EM}$. We note also that this triplet
can have topological implications from the hidden sector. Since $\Delta_{H}$
is an adjoint representation of $SU(2)_{H}$, there exists magnetic monopole
tHooft:1974kcl ; Polyakov:1974ek  and dyon Julia:1975ff  solutions in the
hidden sector which may play the role of topological stable DM Baek:2013dwa .
However, the triplet is not required to generate realistic mass spectra for
all the particles in G2HDM. Therefore if one omits this triplet scalar, the
scalar DM candidate is no longer favorable in the parameter space of G2HDM
according to the analysis in Chen:2019pnt . However, as mentioned above, there
are two other alternative DM candidates in the model. In this paper, we will
show that the non-abelian gauge boson $W^{\prime(p,m)}$ associated with
$SU(2)_{H}$ can be a viable DM as well. In particular we will focus on the low
mass DM scenario in the MeV–GeV range which has attracted a lot of attention
in recent years.

This paper is organized as follows. In Sec. II, we will review some salient
features of the simplified G2HDM without introducing the Higgs triplet field
$\Delta_{H}$ of the extra $SU(2)_{H}$. The theoretical constraints on the
Higgs potential, electroweak precision data for the $Z$-boson mass shift from
Large Electron-Positron Collider (LEP) data, dark photon constraints from
various low energy experiments and the 125 GeV Higgs data constraints from the
Large Hadron Collider (LHC) are discussed in Sec. III. We then turn to the
experimental constraints of dark matter physics in Sec. IV. The cosmological
relic density from Planck satellite Aghanim:2018eyx , underground direct
detection constraints from CRESST III Angloher:2017sxg , DarkSide-50
Agnes:2018ves  and XENON1T Aprile:2019xxb , astrophysical gamma-ray indirect
detection constraints from Fermi-LAT Ackermann:2015zua ; Fermi-LAT:2016uux
and mono-jet constraints from LHC Aaboud:2017phn ; ATLAS:2020wzf ;
Sirunyan:2017hci  for the sub-GeV dark matter $W^{\prime(p,m)}$ in G2HDM are
studied in Secs. IV.1, IV.2, IV.3 and IV.4 respectively. Our numerical results
are presented in Sec. V. We conclude in Sec. VI. In Appendix A, we discuss the
mixing effects in the gauge fixings of the model in a general renormalizable
$R_{\xi}$ gauge. This work can be regarded as an expanded detailed version of
the compact and partial results presented in Ramos:2021omo .

##  II The Simplified G2HDM Model

In this section, we discuss a simplified version of the G2HDM first proposed
in Huang:2015wts . In particular, we will remove the triplet scalar
$\Delta_{H}$ in the original model because it is not absolutely required for a
realistic particle spectra and the number of free parameters in the scalar
potential can be reduced significantly. We note that the Yukawa couplings are
not affected by this simplification since the triplet does not couple to the
fermions in the model.

###  II.1 Particle Content

The gauge group of the simplified G2HDM is the same as in Huang:2015wts ,

${\mathcal{G}}=SU(3)_{C}\times SU(2)_{L}\times U(1)_{Y}\times SU(2)_{H}\times
U(1)_{X}\;.$

In Table 1, we summarize the matter content and their quantum number
assignments in G2HDM. At the minimum risk of confusion, we will continue refer
this model as G2HDM to avoid cluttering throughout the paper with the
adjective word “simplified”. What we are really concerned about is the
electroweak part of $\mathcal{G}$, so the color group $SU(3)_{C}$ is not
relevant in what follows.

Matter Fields | $SU(3)_{C}$ | $SU(2)_{L}$ | $SU(2)_{H}$ | $U(1)_{Y}$ | $U(1)_{X}$ |  $h$-parity
---|---|---|---|---|---|---
$Q_{L}=\left(u_{L}\;\;d_{L}\right)^{\rm T}$ | 3 | 2 | 1 | 1/6 | 0 | $+\;+$
$U_{R}=\left(u_{R}\;\;u^{H}_{R}\right)^{\rm T}$ | 3 | 1 | 2 | 2/3 | 1 | $+\;-$
$D_{R}=\left(d^{H}_{R}\;\;d_{R}\right)^{\rm T}$ | 3 | 1 | 2 | $-1/3$ | $-1$ | $-\;+$
$u_{L}^{H}$ | 3 | 1 | 1 | 2/3 | 0 | $-$
$d_{L}^{H}$ | 3 | 1 | 1 | $-1/3$ | 0 | $-$
$L_{L}=\left(\nu_{L}\;\;e_{L}\right)^{\rm T}$ | 1 | 2 | 1 | $-1/2$ | 0 | $+\;+$
$N_{R}=\left(\nu_{R}\;\;\nu^{H}_{R}\right)^{\rm T}$ | 1 | 1 | 2 | 0 | 1 | $+\;-$
$E_{R}=\left(e^{H}_{R}\;\;e_{R}\right)^{\rm T}$ | 1 | 1 | 2 | $-1$ | $-1$ | $-\;+$
$\nu_{L}^{H}$ | 1 | 1 | 1 | 0 | 0 | $-$
$e_{L}^{H}$ | 1 | 1 | 1 | $-1$ | 0 | $-$
$H=\left(H_{1}\;\;H_{2}\right)^{\rm T}$ | 1 | 2 | 2 | 1/2 | 1 | $+\;-$
$\Phi_{H}=\left(\Phi_{1}\;\;\Phi_{2}\right)^{\rm T}$ | 1 | 1 | 2 | 0 | 1 | $-\;+$
$\mathcal{S}$ | 1 | 1 | 1 | 0 | 0 | $+$
Table 1: Matter content and their quantum number assignments in G2HDM. The
electric charge $Q$ in unit of positron charge $e$ is given by
$Q=T^{3}_{L}+Y$. The scalar $\mathcal{S}$ in the last row is the Stueckelberg
field introduced in Huang:2015wts  to give mass for the $U(1)_{X}$ gauge
boson.

###  II.2 Higgs Potential and Spontaneous Symmetry Breaking

The most general Higgs potential which is invariant under $SU(2)_{L}\times
U(1)_{Y}\times SU(2)_{H}\times U(1)_{X}$ can be written down as follows

$\displaystyle V={}$ $\displaystyle-\mu^{2}_{H}\left(H^{\alpha i}H_{\alpha
i}\right)+\lambda_{H}\left(H^{\alpha i}H_{\alpha
i}\right)^{2}+\frac{1}{2}\lambda^{\prime}_{H}\epsilon_{\alpha\beta}\epsilon^{\gamma\delta}\left(H^{\alpha
i}H_{\gamma i}\right)\left(H^{\beta j}H_{\delta j}\right)$
$\displaystyle-\mu^{2}_{\Phi}\Phi_{H}^{\dagger}\Phi_{H}+\lambda_{\Phi}\left(\Phi_{H}^{\dagger}\Phi_{H}\right)^{2}+\lambda_{H\Phi}\left(H^{\dagger}H\right)\left(\Phi_{H}^{\dagger}\Phi_{H}\right)+\lambda^{\prime}_{H\Phi}\left(H^{\dagger}\Phi_{H}\right)\left(\Phi_{H}^{\dagger}H\right),$
(1)

where ($\alpha$, $\beta$, $\gamma$, $\delta$) and ($i$, $j$) refer to the
$SU(2)_{H}$ and $SU(2)_{L}$ indices respectively, all of which run from 1 to
2, and $H^{\alpha i}=H^{*}_{\alpha i}$.

To facilitate spontaneous symmetry breaking (SSB) and obtain the particle mass
spectra of the model, we shift the fields based on our conventional wisdom

$\displaystyle H_{1}\equiv\begin{pmatrix}H_{11}\\\
H_{12}\end{pmatrix}=\begin{pmatrix}G^{+}\\\
\frac{v+h}{\sqrt{2}}+i\frac{G^{0}}{\sqrt{2}}\end{pmatrix},\;H_{2}\equiv\begin{pmatrix}H_{21}\\\
H_{22}\end{pmatrix}=\begin{pmatrix}H^{+}\\\
H_{2}^{0}\end{pmatrix},\;\Phi_{H}=\begin{pmatrix}G_{H}^{p}\\\
\frac{v_{\Phi}+\phi_{2}}{\sqrt{2}}+i\frac{G_{H}^{0}}{\sqrt{2}}\end{pmatrix},\;$
(2)

where $v$ and $v_{\Phi}$ are the vacuum expectation values (VEV) of $H_{1}$
and $\Phi_{H}$ fields respectively. $H_{2}$ is the inert doublet in G2HDM and
hence does not have VEV. Naively we would think that the Goldstone bosons
$G^{+}$, $G^{p}_{H}$, $G^{0}$ and $G^{0}_{H}$ will be absorbed by the
longitudinal components of $W^{+}$, $W^{\prime p}$, $W^{3}$ and $W^{\prime 3}$
respectively. In Appendix A we will show that the last three Goldstone fields
have mixing effects with other fields in the scalar sector that makes the
situation more interesting but a little bit more complicated.

Substituting the scalar field decomposition of Eq. (2) into the scalar
potential Eq. (II.2) and then minimize the potential, one can obtain the
solutions of VEVs as follows

$\displaystyle v^{2}$ $\displaystyle=$
$\displaystyle\frac{2\left(\lambda_{H\Phi}\mu_{\Phi}^{2}-2\lambda_{\Phi}\mu_{H}^{2}\right)}{\lambda_{H\Phi}^{2}-4\lambda_{H}\lambda_{\Phi}}\;,$
(3) $\displaystyle v_{\Phi}^{2}$ $\displaystyle=$
$\displaystyle\frac{2\left(\lambda_{H\Phi}\mu_{H}^{2}-2\lambda_{H}\mu_{\Phi}^{2}\right)}{\lambda_{H\Phi}^{2}-4\lambda_{H}\lambda_{\Phi}}\;.$
(4)

Equivalently, we can use the minimization conditions to trade $\mu_{H}^{2}$
and $\mu_{\Phi}^{2}$ with $v$ and $v_{\Phi}$ as

$\displaystyle\mu_{H}^{2}$
$\displaystyle=\lambda_{H}v^{2}+\frac{\lambda_{H\Phi}v_{\Phi}^{2}}{2}\;,$ (5)
$\displaystyle\mu_{\Phi}^{2}$
$\displaystyle=\lambda_{\Phi}v_{\Phi}^{2}+\frac{\lambda_{H\Phi}v^{2}}{2}\;.$
(6)

###  II.3 Scalar Mass Spectrum

In the $S=\\{h,\phi_{2}\\}$ basis the mass matrix is given as

${\mathcal{M}}_{S}^{2}=\begin{pmatrix}2\lambda_{H}v^{2}&\lambda_{H\Phi}vv_{\Phi}\\\
\lambda_{H\Phi}vv_{\Phi}&2\lambda_{\Phi}v_{\Phi}^{2}\end{pmatrix}\,.$ (7)

One can use an orthogonal transformation $O^{S}$, which can be parametrized as

$O^{S}=\begin{pmatrix}\cos\theta_{1}&\sin\theta_{1}\\\
-\sin\theta_{1}&\cos\theta_{1}\end{pmatrix}\,,$ (8)

where

$\tan
2\theta_{1}=\frac{2\mathcal{M}^{2}_{S12}}{\mathcal{M}^{2}_{S22}-\mathcal{M}^{2}_{S11}}=\frac{\lambda_{H\Phi}vv_{\Phi}}{\lambda_{\Phi}v_{\Phi}^{2}-\lambda_{H}v^{2}}\;,$
(9)

to diagonalize $\mathcal{M}_{S}^{2}$,

$\left(O^{S}\right)^{\rm T}\cdot{\mathcal{M}}_{S}^{2}\cdot O^{S}={\rm
Diag}\left(m^{2}_{h_{1}},m^{2}_{h_{2}}\right)\;,$ (10)

with $h_{1}$ being identified as the 125 GeV SM-like Higgs boson and $h_{2}$
as a heavier scalar boson. The mass squared eigenvalues of Eq. (7) are given
by

$m_{h_{1,2}}^{2}=\lambda_{H}v^{2}+\lambda_{\Phi}v_{\Phi}^{2}\mp\sqrt{\lambda_{H}^{2}v^{4}-2\lambda_{H}\lambda_{\Phi}v^{2}v_{\Phi}^{2}+\lambda_{H\Phi}^{2}v^{2}v_{\Phi}^{2}+\lambda_{\Phi}^{2}v_{\Phi}^{4}}\;.$
(11)

In the basis of $S^{\prime}=\\{G_{H}^{p},H_{2}^{0*}\\}$, we obtain the mass
matrix:

${\mathcal{M}}_{S^{\prime}}^{2}=\frac{1}{2}\lambda^{\prime}_{H\Phi}\begin{pmatrix}v^{2}&vv_{\Phi}\\\
vv_{\Phi}&v_{\Phi}^{2}\end{pmatrix}\;.$ (12)

Similarly, this mass matrix can be diagonalized by an orthogonal matrix,

$O^{S^{\prime}}=\begin{pmatrix}\cos\theta_{2}&\sin\theta_{2}\\\
-\sin\theta_{2}&\cos\theta_{2}\end{pmatrix}\,,$ (13)

where

$\tan
2\theta_{2}=\frac{2\mathcal{M}^{2}_{S^{\prime}12}}{\mathcal{M}^{2}_{S^{\prime}22}-\mathcal{M}^{2}_{S^{\prime}11}}=\frac{2vv_{\Phi}}{v_{\Phi}^{2}-v^{2}}\;,$
(14)

which gives

$\left(O^{S^{\prime}}\right)^{\rm T}\cdot{\mathcal{M}}_{S^{\prime}}^{2}\cdot
O^{S^{\prime}}={\rm Diag}\left(0,m^{2}_{D}\right)\;.$ (15)

We note that, in Eq. (15), the zero eigenvalue corresponds to the Nambu-
Goldstone boson mass eigenstate $\tilde{G}^{p}_{H}$, while the other
eigenvalue

$m^{2}_{D}=\frac{1}{2}\lambda^{\prime}_{H\Phi}(v^{2}+v_{\Phi}^{2})\;,$ (16)

is the mass of a new dark scalar boson.

The charged Higgs boson mass is given as

$m^{2}_{H^{\pm}}=\frac{1}{2}\left(\lambda^{\prime}_{H\Phi}v_{\Phi}^{2}-\lambda^{\prime}_{H}v^{2}\right)\;.$
(17)

The Goldstone bosons $G^{0}$ and $G^{0}_{H}$ are massless.

The above scalar mass spectrum is derived in the so-called ’t Hooft-Landau
gauge. We will discuss further the mixing effects of the Goldstone bosons with
other scalar fields in a general renormalizable $R_{\xi}$ gauge in Appendix A.

We note that $h_{1,2}$, $G^{0}$ and $G^{0}_{H}$ are even under $h$-parity,
while $\tilde{G}^{p}_{H}$, $D$ and $H^{\pm}$ are odd Chen:2019pnt .

###  II.4 Gauge Sector

After SSB, the $W^{\pm}$ gauge boson of $SU(2)_{L}$ remains the same as in SM
with its mass given by $m_{W}=gv/2$. The $SU(2)_{H}$ gauge boson
$W^{\prime(p,m)}$ receives mass from $\langle H_{1}\rangle$ and
$\langle\Phi_{2}\rangle$ given by

$m_{W^{\prime}}=\frac{1}{2}g_{H}\sqrt{v^{2}+v^{2}_{\Phi}}\;.$ (18)

Note that $W^{\prime(p,m)}$ are electrically neutral and thus do not mix with
the SM $W^{\pm}$. In addition, $W^{\prime(p,m)}$ is odd under $h$-parity. If
it is the lightest $h$-parity odd particle in the model, it will be stable and
can be a DM candidate.

On the other hand, the SM neutral gauge bosons $B$ and $W^{3}$ can mix with
the new gauge bosons $W^{\prime 3}$ and $X$, all of which have even
$h$-parity. Together with the Stueckelberg mass parameters $M_{X}$ and $M_{Y}$
for the two abelian groups $U(1)_{X}$ and $U(1)_{Y}$, SSB generates a $4\times
4$ neutral gauge boson mass matrix in the basis of $\left\\{B,W^{3},W^{\prime
3},X\right\\}~{}$Huang:2015wts ; Huang:2019obt . Due to the theoretical
motivations or prejudices mentioned in Ref. Huang:2019obt , we set the
Stueckelberg mass $M_{Y}=0$. Applying the weak rotation on upper left $2\times
2$ block of the $4\times 4$ mass matrix, one obtains immediately a zero
eigenvalue identified as the SM photon and a $3\times 3$ sub-matrix in the
basis of $\left\\{Z^{\rm SM},W^{\prime 3},X\right\\}$ given by Huang:2019obt ,

${\cal
M}^{2}_{Z}=\begin{pmatrix}m_{Z^{\text{SM}}}^{2}&-\frac{g_{H}v}{2}m_{Z^{\text{SM}}}&-g_{X}vm_{Z^{\text{SM}}}\\\
-\frac{g_{H}v}{2}m_{Z^{\text{SM}}}&m_{W^{\prime}}^{2}&\frac{g_{X}g_{H}\left(v^{2}-v_{\Phi}^{2}\right)}{2}\\\
-g_{X}vm_{Z^{\text{SM}}}&\frac{g_{X}g_{H}\left(v^{2}-v_{\Phi}^{2}\right)}{2}&g_{X}^{2}(v^{2}+v_{\Phi}^{2})+M_{X}^{2}\end{pmatrix}\,,$
(19)

where $g$, $g^{\prime}$, $g_{H}$ and $g_{X}$ are the gauge couplings of
$SU(2)_{L}$, $U(1)_{Y}$, $SU(2)_{H}$ and $U(1)_{X}$ respectively, and
$m_{Z^{\rm SM}}=v\sqrt{g^{2}+g^{\prime 2}}/2$ is the SM $Z$ boson mass
expression. The mass matrix in Eq. (19) can be diagonalized by an orthogonal
rotation matrix ${\cal O}$ so that 111The analytical expression of this
rotation matrix can be found in Ref. Huang:2019obt .

$\begin{pmatrix}Z^{\rm SM}\\\ W^{\prime 3}\\\ X\end{pmatrix}={\cal
O}\cdot\begin{pmatrix}Z\\\ Z^{\prime}\\\ A^{\prime}\end{pmatrix}\;.$ (20)

In this analysis, we arrange the neutral gauge boson masses as
$m_{A^{\prime}}<m_{Z^{\prime}}<m_{Z}\simeq m_{Z^{\rm SM}}$ with $Z$ identified
as the physical $Z$ boson with mass $91.1876\pm 0.0021$ GeV Zyla:2020zbs —the
first heavy neutral vector gauge boson discovered in 1983 at the Super Proton
Synchrotron at CERN. Two more massive neutral vector gauge bosons are
predicted in G2HDM.

By means of a few assumptions motivated by our expectations, it is possible to
draw a few conclusions from Eq. (19). First, the new gauge couplings $g_{H}$
and $g_{X}$ are expected to be much smaller than the SM $g$ and $g^{\prime}$
to avoid large effects on the very precise measurements of the $Z$ properties.
Second, the scale of $v_{\Phi}$ is expected to be larger than $v$ given that
it characterizes the scale of new physics and is directly related to the
masses of beyond the SM (BSM) states. By neglecting any term composed by a
product of any three or more of $g_{H}$, $g_{X}$ and $v^{2}/v_{\Phi}^{2}$ it
is possible to put Eq. (19) into a block diagonal matrix where only $W^{\prime
3}$ and $X$ mix resulting in the approximation

$m_{Z}\approx m_{Z^{\rm SM}}\;.$ (21)

Moreover, the 2$\times$2 squared mass matrix of $W^{\prime 3}$ and $X$ can be
easily diagonalized and somewhat simple approximations can be found. Since we
want the hierarchy $m_{A^{\prime}}<m_{Z^{\prime}}<m_{Z^{\rm SM}}$, the $M_{X}$
parameter is required to have a value smaller than $v$. Assuming $M_{X}<v$
allows us to expand the squared root in the general solution for the
eigenvalues of a 2$\times$2 matrix resulting in the following approximations

$\displaystyle m_{Z^{\prime}}^{2}$ $\displaystyle\approx
m_{W^{\prime}}^{2}\left(1+\frac{4g_{X}^{2}}{g_{H}^{2}}\right)+M_{X}^{2}\left[1-\left(1+\frac{4g_{X}^{2}}{g_{H}^{2}}+\frac{M_{X}^{2}}{m_{W^{\prime}}^{2}}\right)^{-1}\right]\;,$
(22) $\displaystyle m_{A^{\prime}}^{2}$ $\displaystyle\approx
M_{X}^{2}\left(1+\frac{4g_{X}^{2}}{g_{H}^{2}}+\frac{M_{X}^{2}}{m_{W^{\prime}}^{2}}\right)^{-1}\;,$
(23)

where $m_{W^{\prime}}^{2}\approx g_{H}^{2}v_{\Phi}^{2}/4$ was used. From these
expressions we can see that $m_{Z^{\prime}}\gtrsim m_{W^{\prime}}$ and
$m_{A^{\prime}}\lesssim M_{X}$.

Since the couplings of the extra gauge bosons, $Z^{\prime}$ and $A^{\prime}$,
to the SM fermions are proportional to the new gauge couplings $g_{H}$ and/or
$g_{X}$ which are in general much smaller than the SM couplings $g$ and
$g^{\prime}$, the Drell-Yan type processes are suppressed and this can explain
the null results of BSM neutral gauge bosons searches at LEP.

In our study, the couplings of the extra gauge bosons to the SM charged
leptons and to quarks $u$ and $d$ will be important for dark photon
constraints and direct detection of DM. The vectorial and axial parts of their
couplings are given by

$\displaystyle v_{\ell}^{Z(i)}$
$\displaystyle=\mathcal{O}_{1i}\left(-\frac{1}{2}+2s_{W}^{2}\right)-\frac{1}{\sqrt{g^{2}+g^{\prime
2}}}\left(g_{X}\mathcal{O}_{3i}+\frac{1}{2}g_{H}\mathcal{O}_{2i}\right)\;,$
(24) $\displaystyle v_{u}^{Z(i)}$
$\displaystyle=\mathcal{O}_{1i}\left(\frac{1}{2}-\frac{4}{3}s_{W}^{2}\right)+\frac{1}{\sqrt{g^{2}+g^{\prime
2}}}\left(g_{X}\mathcal{O}_{3i}+\frac{1}{2}g_{H}\mathcal{O}_{2i}\right)\;,$
(25) $\displaystyle v_{d}^{Z(i)}$
$\displaystyle=\mathcal{O}_{1i}\left(-\frac{1}{2}+\frac{2}{3}s_{W}^{2}\right)-\frac{1}{\sqrt{g^{2}+g^{\prime
2}}}\left(g_{X}\mathcal{O}_{3i}+\frac{1}{2}g_{H}\mathcal{O}_{2i}\right)\;,$
(26) $\displaystyle a_{\ell}^{Z(i)}$
$\displaystyle=-a_{u}^{Z(i)}=a_{d}^{Z(i)}=-\frac{\mathcal{O}_{1i}}{2}+\frac{1}{\sqrt{g^{2}+g^{\prime
2}}}\left(g_{X}\mathcal{O}_{3i}+\frac{1}{2}g_{H}\mathcal{O}_{2i}\right)\;,$
(27)

where $i=1,2,3$ and $Z(1)\equiv Z$, $Z(2)\equiv Z^{\prime}$, $Z(3)\equiv
A^{\prime}$. Using the relations between mixing matrix elements and mixing
angles in Eqs. (2.6) to (2.9) of Ref. Huang:2019obt , we can find the
following relation

$\frac{2}{\sqrt{g^{2}+g^{\prime
2}}}\frac{\left(g_{X}\mathcal{O}_{3i}+\frac{1}{2}g_{H}\mathcal{O}_{2i}\right)}{O_{1i}}=1-\frac{m_{Z(i)}^{2}}{m_{Z^{\rm
SM}}^{2}}\;,$ (28)

which can be used to simplify the above vectorial and axial couplings to
obtain

$\displaystyle v_{\ell}^{Z(i)}$
$\displaystyle=-\mathcal{O}_{1i}\left(1-2s_{W}^{2}-\frac{r_{i}^{2}}{2}\right)\;,$
(29) $\displaystyle v_{u}^{Z(i)}$
$\displaystyle=\mathcal{O}_{1i}\left(1-\frac{4}{3}s_{W}^{2}-\frac{r_{i}^{2}}{2}\right)\;,$
(30) $\displaystyle v_{d}^{Z(i)}$
$\displaystyle=-\mathcal{O}_{1i}\left(1-\frac{2}{3}s_{W}^{2}-\frac{r_{i}^{2}}{2}\right)\;,$
(31) $\displaystyle a_{\ell}^{Z(i)}$
$\displaystyle=-a_{u}^{Z(i)}=a_{d}^{Z(i)}=-\frac{\mathcal{O}_{1i}}{2}r_{i}^{2}\;,$
(32)

where $r_{i}=m_{Z(i)}/m_{Z^{\rm SM}}$. In this form, it is obvious that the
axial couplings magnitude is expected to be smaller than the ratio of squared
masses $r_{i}^{2}$, which, e.g., for $m_{Z(i)}=1$ GeV is already close to
10-4. In contrast, the vectorial couplings have an $r_{i}$-independent part
whose size is controlled by the mixing matrix element $\mathcal{O}_{1i}$ that
also affects axial couplings. Therefore, for sufficiently light $Z^{\prime}$
and $A^{\prime}$, the contributions from the axial couplings are expected to
be subleading.

The neutral current interactions induced by $W^{\prime(p,m)}$ are given by

${\cal L}\left(W^{\prime}\right)=g_{H}\left(J^{p\,\mu}W_{\mu}^{\prime p}+{\rm
H.c.}\right),$ (33)

where

$J^{p\,\mu}=\frac{1}{\sqrt{2}}\left[\overline{u_{R}}V_{u}^{H}\gamma^{\mu}u_{R}^{H}+\overline{d_{R}^{H}}\left(V_{d}^{H}\right)^{\dagger}\gamma^{\mu}d_{R}+\overline{\nu_{R}}V_{\nu}^{H}\gamma^{\mu}\nu_{R}^{H}+\overline{e_{R}^{H}}\left(V_{e}^{H}\right)^{\dagger}\gamma^{\mu}e_{R}\right]\;,$
(34)

with $V_{u}^{H}$, $V_{d}^{H}$, $V_{\nu}^{H}$ and $V_{e}^{H}$ being the new
unitary mass rotation matrices for the fermions. The neutral current
interactions induced by $Z$, $Z^{\prime}$, $A^{\prime}$ bosons can be found in
Ref. Huang:2019obt .

###  II.5 Free Parameters

From the scalar potential we recognize 9 free parameters including couplings
and VEVs. Of those 9 parameters, $\mu_{H}^{2}$ and $\mu_{\Phi}^{2}$ can be
related to other parameters using the minimization conditions according to
Eqs. (5) and (6), leaving only 7 free parameters. Furthermore, Eq. (10) can be
used to relate the parameters $\lambda_{H}$, $\lambda_{\Phi}$ and
$\lambda_{H\Phi}$ to the physical squared masses $m_{h_{1}}^{2}$,
$m_{h_{2}}^{2}$ and the mixing angle $\theta_{1}$ obtaining the following
relations

$\displaystyle\lambda_{H}$
$\displaystyle=\frac{1}{2v^{2}}\left(m^{2}_{h_{1}}\cos^{2}{\theta_{1}}+m^{2}_{h_{2}}\sin^{2}{\theta_{1}}\right)\;,$
(35) $\displaystyle\lambda_{\Phi}$
$\displaystyle=\frac{1}{2v_{\Phi}^{2}}\left(m^{2}_{h_{1}}\sin^{2}{\theta_{1}}+m^{2}_{h_{2}}\cos^{2}{\theta_{1}}\right)\;,$
(36) $\displaystyle\lambda_{H\Phi}$
$\displaystyle=\frac{1}{2vv_{\Phi}}\left[\left(m^{2}_{h_{2}}-m^{2}_{h_{1}}\right)\sin{\left(2\theta_{1}\right)}\right]\;.$
(37)

Using the mass of $D$, Eq. (16), we can express $\lambda^{\prime}_{H\Phi}$ in
the form

$\lambda^{\prime}_{H\Phi}=\frac{2m^{2}_{D}}{v^{2}+v_{\Phi}^{2}}\;.$ (38)

Using this last expression and the mass of the charged Higgs of Eq. (17), we
can write $\lambda^{\prime}_{H}$ as

$\lambda^{\prime}_{H}=\frac{2}{v^{2}}\left[\frac{m^{2}_{D}v_{\Phi}^{2}}{v^{2}+v_{\Phi}^{2}}-m^{2}_{H^{\pm}}\right]\;.$
(39)

Finally, we can use Eq. (18) to relate $v_{\Phi}$ to the $W^{\prime}$ mass as

$v_{\Phi}^{2}=\frac{4m_{W^{\prime}}^{2}}{g_{H}^{2}}-v^{2}\;.$ (40)

Using the expressions in this subsection allow us to trade six model
parameters with five physical squared masses and one mixing angle,

$\lambda_{H},\lambda_{\Phi},\lambda_{H\Phi},\lambda^{\prime}_{H\Phi},\lambda^{\prime}_{H},v_{\Phi}\>\>\>\>\rightarrow\>\>\>\>m^{2}_{h_{1}},m^{2}_{h_{2}},m^{2}_{D},m^{2}_{H^{\pm}},m^{2}_{W^{\prime}},\theta_{1}\;.$
(41)

The remaining free parameters of the model are the heavy fermion masses
$m_{f^{\rm H}}$, the Stueckelberg mass $M_{X}$, and the gauge couplings
$g_{X}$ and $g_{H}$. Considering that the mass of the Higgs, $m_{h_{1}}$, has
already been measured Zyla:2020zbs , we are left with a total of 8 free
parameters plus the masses for 12 heavy hidden fermions. The effects of heavy
hidden fermions in complex scalar dark matter phenomenology in G2HDM was
analyzed in Dirgantara:2020lqy .

##  III Constraints

In this section, we examine the model using various constraints including the
theoretical constraints of the scalar potential, electroweak precision data,
dark photon physics and Higgs measurements at LHC. DM constraints are
presented separately in the next section.

###  III.1 Theoretical Constraints

The theoretical constraints on the original G2HDM were studied in Ref.
Arhrib:2018sbz . Here, we follow closely on the steps of that work but remove
the scalar triplet and related parameters from the original model.

Vacuum Stability: To make sure the scalar potential is bounded from below, the
sum of all quartic terms in the scalar potential needs to be positive. In the
same way as Ref. Arhrib:2018sbz  we use copositivity conditions given by the
following constraints

$\widetilde{\lambda}_{H}(\eta)\geq 0\;,\;\;\;\lambda_{\Phi}\geq
0\;,\;\;\;\widetilde{\lambda}_{H\Phi}(\xi)+2\sqrt{\widetilde{\lambda}_{H}(\eta)\lambda_{\Phi}}\geq
0\;,$ (42)

where
$\widetilde{\lambda}_{H}(\eta)\equiv\lambda_{H}+\eta\lambda^{\prime}_{H}$ and
${\widetilde{\lambda}_{H\Phi}}(\xi)\equiv\lambda_{H\Phi}+\xi\lambda^{\prime}_{H\Phi}$.
The conditions of Eq. (42) have to be met for any value of $\xi$ and $\eta$ in
the ranges $0\leq\xi\leq 1$ and $-1\leq\eta\leq 0$.

Perturbative Unitarity: It is required that the parameter space remains within
the perturbative limits. To this end, here we compute the 2$\rightarrow$2
scalar scattering amplitudes induced by the quartic couplings. The
2$\rightarrow$2 processes induced by vertices from scalar cubic couplings and
gauge interactions are suppressed by large momentum exchange in their
propagators Arhrib:2018sbz .

Perturbative unitarity requires

$\displaystyle|\lambda_{H}|,|\lambda_{\Phi}|\leq 4\pi\;,|\lambda_{H\Phi}|\leq
8\pi\;,|\lambda^{\prime}_{H\Phi}|\;,|\lambda_{H}^{\prime}|\leq
8\sqrt{2}\pi\;,$ (43) $\displaystyle|2\lambda_{H}\pm\lambda^{\prime}_{H}|\leq
8\pi\;,|\lambda_{H\Phi}+\lambda_{H\Phi}^{\prime}|\leq 8\pi\;,$ (44)
$\displaystyle\left|(\lambda_{H}+\lambda_{H}^{\prime}/2+\lambda_{\Phi})\pm\sqrt{2\lambda^{\prime
2}_{H\Phi}+(\lambda_{H}+\lambda_{H}^{\prime}/2-\lambda_{\Phi})^{2}}\right|\leq
8\pi\;,$ (45)
$\displaystyle\left|(5\lambda_{H}-\lambda_{H}^{\prime}/2+3\lambda_{\Phi})\pm\sqrt{(5\lambda_{H}-\lambda_{H}^{\prime}/2-3\lambda_{\Phi})^{2}+2(2\lambda_{H\Phi}+{\lambda}^{\prime}_{H\Phi})^{2}}\right|\leq
8\pi\;.$ (46)

###  III.2 Electroweak Constraints

A comprehensive study on the electroweak precision constraints in the original
model has been preformed in Ref. Huang:2019obt . We found that the $Z$ mass
shift is the most stringent among all the electroweak precision constraints.
In particular, for the parameter space of interest in our analysis, _i.e._ ,
$g_{H}\sim g_{X}\ll 1$ and $(m_{W^{\prime}}^{2},M_{X}^{2})\ll m_{Z^{\rm
SM}}^{2}$, the $Z$ mass shift, $|\Delta m_{Z}|=|m_{Z}-m_{Z^{\rm SM}}|$, can be
estimated as follows

$\left|\frac{\Delta m_{Z}}{m_{Z^{\rm
SM}}}\right|\simeq\frac{5}{2}\frac{g_{H}^{2}}{g_{Z}^{2}}\left(1+\frac{7}{5}\frac{m_{W^{\prime}}^{2}}{m_{Z^{\rm
SM}}^{2}}-\frac{4}{5}\frac{M_{X}^{2}}{m_{Z^{\rm SM}}^{2}}\right)^{-1}\;,$ (47)

where $g_{Z}=g/\cos\theta_{W}$ with $\theta_{W}$ being the Weinberg angle.
Following the methodology of Ref. Feldman:2006ce , we can obtain the
experimental uncertainty of the $Z$ mass as

$\bigg{[}\frac{\delta m_{Z}}{m_{Z^{\rm
SM}}}\bigg{]}^{2}=\bigg{[}\frac{c_{W}^{-2}-2t_{W}^{2}}{\delta
m_{W}^{-1}m_{W}}\bigg{]}^{2}+\;\frac{t_{W}^{4}(\delta\Delta r)^{2}}{4(1-\Delta
r)^{2}}\;,$ (48)

where $t_{W}=\tan\theta_{W}$ and $\Delta r$ is the radiative correction. Using
the PDG values Zyla:2020zbs  of $m_{W}\pm\delta m_{W}=80.387\pm 0.016$ GeV,
$\Delta r\pm\delta\Delta r=0.03652\mp 0.00021\pm 0.00007$ and
$\sin\theta_{W}=0.22343$, and by requiring $|\Delta m_{Z}|<|\delta m_{Z}|$,
one obtains an upper bound on $g_{H}$ and $g_{X}$

$|g_{X}|\sim|g_{H}|\lesssim
0.006\,\times\sqrt{1-\frac{7}{5}\frac{m_{W^{\prime}}^{2}}{m_{Z^{\rm
SM}}^{2}}+\frac{4}{5}\frac{M_{X}^{2}}{m_{Z^{\rm SM}}^{2}}}\;.$ (49)

###  III.3 Dark Photon

The light boson $A^{\prime}$ can be treated as a dark photon and, therefore,
dark photon constraints have to be applied. In particular, due to the vertex
$A^{\prime}\bar{\ell}\ell$, where $\ell$ represents any charged lepton, it is
expected that for a sufficiently large coupling it should be possible to
observe the $A^{\prime}$ resonance in the invariant mass distribution of
$e^{+}e^{-}$ and $\mu^{+}\mu^{-}$. Dark photon experiments constrain the size
of the coupling via a parameter $\varepsilon_{\ell}$. In the decay width
$\Gamma\left(A^{\prime}\rightarrow\ell\ell\right)$ the parameter
$\varepsilon_{\ell}$ appears as Fabbrichesi:2020wbt

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

where $\mu_{\ell}=2m_{\ell}/m_{A^{\prime}}<1$ since this decay channel only
opens for $m_{A^{\prime}}>2m_{\ell}$. In the G2HDM, the parameter
$\varepsilon_{\ell}$ at tree level is given by

$\varepsilon_{\ell}=\frac{1}{2s_{W}c_{W}}\sqrt{\left(v^{A^{\prime}}_{\ell}\right)^{2}+\left(a^{A^{\prime}}_{\ell}\right)^{2}\left(\frac{1-\mu_{\ell}^{2}}{1+\mu_{\ell}^{2}/2}\right)}\;,$
(51)

where $v^{A^{\prime}}_{\ell}$ and $a^{A^{\prime}}_{\ell}$ are given in Sec. 2
of Ref. Huang:2019obt . From Ref. Huang:2019obt  we know that
$v_{\ell}^{A^{\prime}}$ and $a_{\ell}^{A^{\prime}}$ are the same for all the
charged leptons and, thus, the only distinction in $\varepsilon_{\ell}$
between different flavors of leptons comes from $\mu_{\ell}$. As mentioned in
Sec. II.4, for light enough $A^{\prime}$ the axial coupling will be negligible
and $\varepsilon_{\ell}$ is expected to be nearly independent of $\mu_{\ell}$,
as is usually the case in models with dark photon. It is important to mention
that, since $Z^{\prime}$ is also expected to be light, the dark photon
experimental limits can also be applied as above with $A^{\prime}\rightarrow
Z^{\prime}$ in Eq. (50) and
$\\{v_{\ell}^{A^{\prime}},a_{\ell}^{A^{\prime}}\\}\rightarrow\\{v^{Z^{\prime}}_{\ell},a^{Z^{\prime}}_{\ell}\\}$
in Eq. (51). However, since $A^{\prime}$ is lighter by definition it is
expected to be more strongly constrained.

There are several experiments with reported stringent limits for
$m_{A^{\prime}}>1$ MeV Aaij:2019bvg ; Lees:2014xha ; Batley:2015lha ;
Banerjee:2018vgk ; Riordan:1987aw ; Blumlein:2011mv ; Blumlein:2013cua ;
darkphotonexp . The existing limits on $\varepsilon_{\ell}$ for a dark photon
mass $m_{A^{\prime}}>1$ MeV are displayed on the top pane of Fig. 10 in Ref.
Fabbrichesi:2020wbt .

###  III.4 Higgs Collider Data

####  III.4.1 Higgs Boson Mass

As aforementioned, $h_{1}$ is identified as the observed Higgs boson at the
LHC. In this analysis, we take the mass of Higgs boson as $m_{h_{1}}=125.10\pm
0.14$ GeV Zyla:2020zbs .

####  III.4.2 Higgs Decays into Diphoton

The decay rate for $h_{1}\rightarrow\gamma\gamma$ is given by

$\Gamma(h_{1}\rightarrow\gamma\gamma)=\frac{1}{64\pi}m_{h_{1}}^{3}\biggl{|}F_{\gamma\gamma}(W^{\pm})+F_{\gamma\gamma}(H^{\pm})+\sum_{{\rm
Charged}\,f^{\rm SM}}F_{\gamma\gamma}(f^{\rm SM})+\sum_{{\rm
Charged}\,f^{H}}F_{\gamma\gamma}(f^{\rm H})\biggr{|}^{2}\;,$ (52)

where 222 Note that both the charged Higgs and heavy fermion contributions in
$h_{1}\rightarrow\gamma\gamma$ were not handled properly in Huang:2015wts .

$\displaystyle F_{\gamma\gamma}(W^{\pm})={}$
$\displaystyle\frac{-1}{16\pi^{2}}\cdot e^{2}\cdot
g\cdot\frac{1}{m_{W}}\cdot\cos\theta_{1}\cdot\left[2+3\tau_{W}+3\tau_{W}\left(2-\tau_{W}\right)f(\tau_{W})\right]\,,$
(53) $\displaystyle F_{\gamma\gamma}(H^{\pm})={}$
$\displaystyle\frac{-1}{16\pi^{2}}\cdot e^{2}\cdot
g_{h_{1}H^{+}H^{-}}\cdot\frac{1}{m^{2}_{H^{\pm}}}\cdot\left\\{\tau_{H^{\pm}}\left[1-\tau_{H^{\pm}}f(\tau_{H^{\pm}})\right]\right\\}\;,$
(54) $\displaystyle F_{\gamma\gamma}(f^{\rm SM})={}$
$\displaystyle\frac{1}{16\pi^{2}}\cdot N_{c}\cdot e^{2}Q^{2}_{f^{\rm
SM}}\cdot\frac{4}{v}\cdot\cos\theta_{1}\cdot\left\\{\tau_{f^{\rm
SM}}\left[1+\left(1-\tau_{f^{\rm SM}}\right)f(\tau_{f^{\rm
SM}})\right]\right\\}\;,$ (55) $\displaystyle F_{\gamma\gamma}(f^{\rm H})={}$
$\displaystyle\frac{-1}{16\pi^{2}}\cdot N_{c}\cdot e^{2}Q^{2}_{f^{\rm
H}}\cdot\frac{4}{v_{\Phi}}\sin\theta_{1}\cdot\left\\{\tau_{f^{\rm
H}}\left[1+\left(1-\tau_{f^{\rm H}}\right)f(\tau_{f^{\rm
H}})\right]\right\\}\;.$ (56)

Here, and in what follows, we define $\tau_{i}=4\,m_{i}^{2}/m_{h_{1}}^{2}$
where $i$ indicates which particle is running inside the loop. $N_{c}$ is the
color factor, 1 for leptons and 3 for quarks. The coupling
$g_{h_{1}H^{+}H^{-}}$ corresponds to the $h_{1}H^{+}H^{-}$ vertex and is given
by

$g_{h_{1}H^{+}H^{-}}=\left(2\lambda_{H}-\lambda^{\prime}_{H}\right)v\cos\theta_{1}-\left(\lambda_{H\Phi}+\lambda^{\prime}_{H\Phi}\right)v_{\Phi}\sin\theta_{1}\;.$
(57)

The well-known loop function $f(x)$ is

$f(x)=\left\\{\begin{array}[]{lr}\arcsin^{2}\left(\frac{1}{\sqrt{x}}\right)&(x\geq
1)\,,\\\
-\frac{1}{4}\left[\ln\left(\frac{1+\sqrt{1-x}}{1-\sqrt{1-x}}\right)-i\pi\right]^{2}&(x<1)\,.\end{array}\right.$
(58)

The signal strength parameter for the Higgs boson produced from the gluon-
gluon fusion (ggH) can be obtained as

$\mu^{\gamma\gamma}_{\rm ggH}=\frac{\Gamma_{h}^{\rm
SM}}{\Gamma_{h_{1}}}\frac{\Gamma({h_{1}\rightarrow
gg})\Gamma({h_{1}\rightarrow\gamma\gamma})}{\Gamma^{\rm SM}({h\rightarrow
gg})\Gamma^{\rm SM}({h\rightarrow\gamma\gamma})}\;,$ (59)

where the superscript SM refers to the SM Higgs boson $h$. The decay width of
$h_{1}$ into two gluons is given by Huang:2015wts

$\Gamma(h_{1}\rightarrow
gg)=\frac{\alpha_{s}^{2}}{48\pi^{3}}m_{h_{1}}^{3}\left|\sum_{q^{\rm
SM}}F_{gg}(q^{\rm SM})+\sum_{q^{\rm H}}F_{gg}(q^{\rm H})\right|^{2}\;,$ (60)

where

$\displaystyle F_{gg}(q^{\rm SM})={}$
$\displaystyle\frac{\cos\theta_{1}}{v}\tau_{q^{\rm
SM}}\left[1+\left(1-\tau_{q^{\rm SM}}\right)f(\tau_{q^{\rm SM}})\right],$ (61)
$\displaystyle F_{gg}(q^{\rm H})={}$
$\displaystyle-\frac{\sin\theta_{1}}{v_{\Phi}}\tau_{q^{\rm
H}}\left[1+\left(1-\tau_{q^{\rm H}}\right)f(\tau_{q^{\rm H}})\right]\;,$ (62)

with $q^{\rm SM}$ and $q^{\rm H}$ refer to the SM quarks and the new colored
fermions. The latest measurement of this signal strength is given by ATLAS as
0.96$\pm$0.14 Aad:2019mbh .

####  III.4.3 $h_{1}\rightarrow$ SM Fermions

The decay width of Higgs boson to SM fermions is given by

$\Gamma({h_{1}\rightarrow f^{\rm SM}\bar{f}^{\rm
SM}})=\frac{N_{c}}{8\pi}\frac{m_{h_{1}}m_{f^{\rm
SM}}^{2}}{v^{2}}\left(1-\tau_{f^{\rm SM}}\right)^{3/2}\cos^{2}\theta_{1}\;.$
(63)

The signal strength for ggH production is then given by

$\mu^{ff}_{\rm ggH}=\cos^{2}\theta_{1}\frac{\Gamma(h^{\rm
SM})}{\Gamma_{h_{1}}}\frac{\Gamma({h_{1}\rightarrow gg})}{\Gamma^{\rm
SM}({h\rightarrow gg})}\;.$ (64)

Note that this expression is independent of the fermion flavor and thus we
compare against the best measured signal strength given by the decay into a
pair of $\tau^{+}\tau^{-}$: $\mu^{\tau\tau}_{\rm ggH}=1.05^{+0.53}_{-0.47}$
Sirunyan:2018koj .

####  III.4.4 Invisible Higgs Decay

If $m_{h_{1}}>2\,m_{W^{\prime}}$, the Higgs boson can decay invisibly into a
pair of $W^{\prime(p,m)}$. The decay width of $h_{1}\rightarrow W^{\prime
p}W^{\prime m}$ is given by

$\Gamma({h_{1}\rightarrow W^{\prime p}W^{\prime
m}})=\frac{g_{H}^{4}\left(v\cos\theta_{1}-v_{\Phi}\sin\theta_{1}\right)^{2}}{256\pi}\frac{m_{h_{1}}^{3}}{m_{W^{\prime}}^{4}}\left(1-\tau_{W^{\prime}}+\frac{3}{4}\tau_{W^{\prime}}^{2}\right)\sqrt{1-\tau_{W^{\prime}}}\;.$
(65)

In our parameter choice, we will assume $2\,m_{\nu^{H}}>m_{h_{1}}$ so that
$h_{1}$ does not decay into a pair of $\nu^{H}$. The branching ratio of
invisible Higgs decay is then given by

${\rm BR}({h_{1}\rightarrow{\rm inv}})=\frac{\Gamma({h_{1}\rightarrow
W^{\prime p}W^{\prime m}})}{\Gamma_{h_{1}}}\;.$ (66)

Recently, the ATLAS collaboration reported the most stringent constraint on
the invisible decays of the Higgs produced via vector boson fusion. Assuming
that the Higgs boson production cross section is comparable to the SM, the
ATLAS collaboration set the limit ${\rm BR}({h_{1}\rightarrow{\rm inv}})<0.13$
at $95\%$ C.L. ATLAS:2020cjb .

##  IV Dark Matter Constraints

###  IV.1 Relic Density

The DM scenario presented here works similarly to the very well known WIMP DM.
The DM candidate $W^{\prime(p,m)}$ begins in thermal equilibrium with the
particle species in the hot primordial soup in the early universe before
starting to freeze-out due to the expansion of the universe. The Boltzmann
equation allows us to determine the evolution of the DM density and precisely
determine the amount of relics that remain after freeze-out. This evolution is
heavily influenced by the back and forth annihilation (creation) of pair of DM
particles into (from) pairs of SM states in the early universe and their
number densities. An excessive annihilation of DM into SM particles would
result in very low relic density while not having enough annihilation would
leave an overabundant DM. In our model, couplings between SM and BSM states
have to be suppressed to minimize the effects on the precisely measured
properties of the $Z$ and the Higgs, given that these measurements are in good
agreement with the SM.

Figure 1:  Dominant annihilation channels of the $W^{\prime\,(p,m)}$ DM
candidate.

The main DM annihilation channels in our model are to pairs of SM fermions
mediated by $Z$ and $Z^{\prime}$ as depicted in Fig. 1. Other annihilation
channels are also possible but are far more suppressed compared to the
channels just mentioned. First, there is the $A^{\prime}$ exchange diagram.
The $A^{\prime}$ couplings to SM fermions are suppressed by combinations of
new gauge couplings in $v^{A^{\prime}}_{f}$ and $a^{A^{\prime}}_{f}$. Similar
to the case of $A^{\prime}$, the $Z^{\prime}$ couplings to the SM fermions are
suppressed by its own $v^{Z^{\prime}}_{f}$ and $a^{Z^{\prime}}_{f}$. However,
it is possible to have $m_{Z^{\prime}}$ close to twice the mass of the
$W^{\prime(p,m)}$, resulting in an important contribution from resonant
annihilations. Secondly, we also have the $h_{1}$ and $h_{2}$ Higgs exchange
diagrams. Their couplings to pairs of $W^{\prime(p,m)}$ and SM fermions are
suppressed by $g_{H}$ and light fermion masses $m_{q}/v$ respectively.
Finally, it is possible to have $t$-channel annihilation diagram via the
exchange of a new heavy fermion, $f^{\rm H}$, but this channel is suppressed
by a factor of $g_{H}$ on each of the two vertices of the diagram and by the
mass of the heavy fermion in the propagator. The Feynman diagrams for the main
annihilation processes as given in Fig. 1 can be computed straightforwardly
and give rise to the following total cross section for each final fermion pair

$\displaystyle\sigma(W^{\prime p}W^{\prime m}\rightarrow\bar{f}f)={}$
$\displaystyle\frac{N_{c}g_{M}^{2}g_{H}^{2}m_{W^{\prime}}^{2}}{72\pi
s^{2}}\left(1-\frac{4\,m_{W^{\prime}}^{2}}{s}\right)^{1/2}\left(1-\frac{4\,m_{f}^{2}}{s}\right)^{1/2}\left(\frac{s^{2}}{m_{W^{\prime}}^{4}}+20\frac{s}{m_{W^{\prime}}^{2}}+12\right)$
$\displaystyle\times\left(\frac{\left(s-m_{f}^{2}\right)}{6m_{W^{\prime}}^{2}}\mathcal{V}_{+}+\frac{m_{f}^{2}}{2m_{W^{\prime}}^{2}}\mathcal{V}_{-}\right)\,,$
(67)

where

$\displaystyle g_{M}=\frac{\sqrt{g^{2}+g^{\prime
2}}}{2}\,,\qquad\mathcal{\hat{D}}_{k}=1-\frac{m_{k}^{2}}{s}+i\frac{\Gamma_{k}m_{k}}{s}\,,$
(68) $\displaystyle\mathcal{V}_{\pm}=\mathcal{O}_{21}^{2}\frac{v_{f}^{2}\pm
a_{f}^{2}}{|\hat{\mathcal{D}}_{Z}|^{2}}+\mathcal{O}_{22}^{2}\frac{(v_{f}^{Z^{\prime}})^{2}\pm(a_{f}^{Z^{\prime}})^{2}}{|\hat{\mathcal{D}}_{Z^{\prime}}|^{2}}+2\mathcal{O}_{22}\mathcal{O}_{21}\left(v_{f}v^{Z^{\prime}}_{f}\pm
a_{f}a^{Z^{\prime}}_{f}\right){\rm
Re}\\!\left(\frac{1}{\hat{\mathcal{D}}_{Z}\hat{\mathcal{D}}^{*}_{Z^{\prime}}}\right)\;.$
(69)

The annihilation mediated by the $Z$ is suppressed by a factor of
${\mathcal{O}}_{21}^{2}$ required to be small mostly by measurements on the
decay width of the $Z$ and the decay branching fractions that limit the
$Z\rightarrow W^{\prime p}W^{\prime m}$ process. As mentioned above, the
channel mediated by the $Z^{\prime}$ is also suppressed by combinations of
gauge couplings that make most of the size of $v^{Z^{\prime}}_{f}$ and
$a^{Z^{\prime}}_{f}$. However, these suppressions are not as strong as the
suppression in other channels and when we include the effects from
$Z^{\prime}$ resonance it is possible to bring the relic density to its
expected value of $\Omega h^{2}\sim 0.1$.

In our study, we will consider the measured value of $\Omega h^{2}=0.120\pm
0.001$ as given by the Planck collaboration Aghanim:2018eyx .

###  IV.2 Direct Detection

Due to the small coupling between the DM candidate, $W^{\prime(p,m)}$, with
the SM-like states $h_{1}$ and $Z$ and the BSM states $h_{2}$, $Z^{\prime}$
and $A^{\prime}$ which couple to the visible sector, it is possible to have
effects from DM scattering against nucleons in detectors used in direct
detection experiments. In this case, we have to consider the elastic
scattering between a DM particle and the partons (both quarks and gluons)
present in the nucleon. The suppression of vertices works in the same way as
in the annihilation processes described in Sec. IV.1, where $h_{j}$ mediated
processes are suppressed by a factor of $g_{H}^{2}m_{q}^{2}/v^{2}$ in the
cross section with an additional $m_{h_{j}}^{-4}$ suppression from the
propagator since these interactions happen via $t$-channel. Therefore, we are
only left with the processes mediated by $Z$, $Z^{\prime}$ and $A^{\prime}$ in
the $t$-channel. Usually, for direct detection processes the momentum exchange
is considered to be very small and therefore $t$ is expected to be small as
well. This will result in amplitudes suppressed by the inverse squared of the
mass of the mediator meaning that the light states, $Z^{\prime}$ and
$A^{\prime}$, will be less suppressed. In the approximation where the momentum
exchange is smaller than the mass of the mediator, we can write the
interaction between DM and light quark $q$ as a contact interaction given by

$\mathcal{L}_{\rm CI-
DD}=\sum_{q}\sum_{i=2}^{3}\frac{g_{M}g_{H}\mathcal{O}_{2i}v^{Z(i)}_{q}}{2m_{Z(i)}^{2}}\left(W^{\prime
p\,\mu}\partial_{\nu}W^{\prime m}_{\mu}-W^{\prime
m\,\mu}\partial_{\nu}W^{\prime p}_{\mu}\right)\bar{q}\gamma^{\nu}q\;,$ (70)

where $i=\\{2,3\\}$ corresponds to $Z(2)\equiv Z^{\prime}$ and $Z(3)\equiv
A^{\prime}$ respectively. It is worth noting that, as light as the mediators
$Z^{\prime}$ and $A^{\prime}$ are, we can still integrate them out thanks to
the comparably small maximum momentum transfer, $q_{\max}$. The smallness of
$q_{\max}$ is mostly due to $W^{\prime}$ being small as well. Consider
$q_{\max}\sim 2\,v_{\rm DM}\,m_{W^{\prime}}m_{A}/(m_{W^{\prime}}+m_{A})$ with
$m_{W^{\prime}}=0.5$ GeV and $v_{\rm DM}=10^{-3}$ c, and the target mass
$m_{A}=131$ GeV or $40$ GeV for xenon or argon target respectively. In both
cases $q_{\max}\sim O(1\ \text{MeV})$ while we expect $m_{A^{\prime}}\gtrsim
O(10\ \text{MeV})$ due to constraints on dark photons. Additionally, smaller
$m_{W^{\prime}}$ results in even smaller $q_{\max}$. Furthermore, for the
axial part of the interaction with the quark, in the small momentum exchange
limit, only the space components of $\gamma^{\nu}$ remain but these components
are suppressed by the $W^{\prime(p,m)}$ momentum due to the derivatives
$\partial_{\nu}W^{\prime(p,m)}$ in Eq. (70) Arcadi:2017kky ; Escudero:2016gzx
. This, together with axial couplings that are comparably much smaller than
the vectorial ones results in an spin dependent cross section that is expected
to be several orders of magnitude smaller than the spin independent one.

From Eq. (70), it is clear that the $A^{\prime}$ mediated process is expected
to dominate the cross section unless
$|\mathcal{O}_{23}/\mathcal{O}_{22}|<|m_{A^{\prime}}/m_{Z^{\prime}}|^{2}$. The
case where both mediators participate equally is expected to happen only
through fine tuning of masses and mixings. Therefore, we expect the cross
section with the nucleons to be mostly mediated by either $A^{\prime}$ or
$Z^{\prime}$. The elastic cross section between $W^{\prime(p,m)}$ and a
nucleon, $N$, is given by

$\displaystyle\sigma^{\rm SI}_{W^{\prime}N}$ $\displaystyle=\sigma^{\rm
SI}_{W^{\prime}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}}\;,$
(71) $\displaystyle\sigma^{\rm SI}_{W^{\prime}p}$
$\displaystyle=\frac{\mu_{p}^{2}g_{M}^{2}g_{H}^{2}\mathcal{O}_{2i}^{2}}{4\pi
m_{Z(i)}^{4}}f_{p}^{2}\;,$ (72)

where $i=2$ or $3$ depending on the dominant mediator according to the
discussion above, $\mu_{p}=m_{W^{\prime}}m_{p}/(m_{W^{\prime}}+m_{p})$ is the
reduced DM-proton mass,
$\mu_{A_{k}}=m_{W^{\prime}}m_{A_{k}}/(m_{W^{\prime}}+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 $Z_{\text{atom}}$ and the isotope dependent variables $\eta_{k}$ and
$A_{k}$ are the abundance and mass number of the $k^{\text{th}}$ target
isotope, respectively. Direct detection experiments usually report the number
in Eq. (71) assuming isospin conservation, _i.e._ , $f_{p}=f_{n}$. In that
case, it is straightforward to see that the ratio of the sums over isotopes
reduces to 1 and $\sigma^{\rm SI}_{W^{\prime}N}=\sigma^{\rm
SI}_{W^{\prime}p}$. However, in our case the couplings between quarks, $u$ and
$d$, and the gauge bosons, $Z^{\prime}$ and $A^{\prime}$, are all different
due to their distinct SM charges leading to isospin violation (ISV), _i.e._ ,
$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}_{W^{\prime}p}/\sigma^{\rm SI}_{W^{\prime}N}$ to account for ISV effects
and use it to limit $\sigma^{\rm SI}_{W^{\prime}p}$ as given by Eq. (72). This
rescaling depends on the mass of DM, the atomic numbers and the ratio
$f_{n}/f_{p}$, and, therefore, will be different for different points in the
parameter space.

To constraint the $W^{\prime(p,m)}$-proton cross section we will use the most
recent upper limits set by the experiments CRESST III Angloher:2017sxg ,
DarkSide-50 Agnes:2018ves  and XENON1T Aprile:2019xxb .

###  IV.3 Indirect Detection

Due to DM annihilation before freeze out happening through the resonance of an
otherwise suppressed channel, the annihilation of DM in the present—after the
shift in energy from the early to the current Universe—loses the resonance
resulting in a very low annihilation cross section. We have checked that the
value of the total annihilation cross section in G2HDM at the present time is
of order $10^{-32}$ cm${}^{3}\cdot$s-1 or below, much lower than the canonical
limits set for various channels by Fermi-LAT data Ackermann:2015zua ; Fermi-
LAT:2016uux .

###  IV.4 Mono-Jet

BP | $m_{W^{\prime}}$ | $M_{X}$ | $m_{h_{2}}$ | $m_{D}$ | $m_{H^{\pm}}$ | $g_{H}$ | $g_{X}$ | $\theta_{1}$ | $\sigma_{pp\rightarrow W^{\prime p}W^{\prime m}j}^{\rm precut}$
---|---|---|---|---|---|---|---|---|---
| (GeV) | (GeV) | (TeV) | (TeV) | (TeV) | ($10^{-4}$) | ($10^{-4}$) | (rad) | (fb)
1 | 1.0 | 1.67 | 0.79 | 2.0 | 2.12 | 5.4 | 3.2 | 0.17 | 3.0
2 | 0.17 | 0.33 | 2.9 | 1.54 | 1.64 | 1.0 | 0.4 | 0.11 | 3.8
Table 2: Parameters for the two benchmark points where the production cross
section of the mono-jet signals with precuts are computed, as shown in the
last column.

The occurrence of energetic jets with large missing transverse momentum has
been searched by ATLAS Aaboud:2017phn ; ATLAS:2020wzf  and CMS
Sirunyan:2017hci  collaborations. However, the observed results are overall in
agreement with the SM predictions and only exclusion limits have been
reported. In our model, the process $pp\rightarrow W^{\prime p}W^{\prime m}j$
can give rise to mono-jet events at the LHC.

To analyze the mono-jet signal at the LHC, we choose two benchmark points
(BPs) which are shown in Table 2. These two BPs satisfy all theoretical, Higgs
data and DM constraints. We evaluate the signal process cross section using
MadGraph 5 Alwall:2011uj  with precuts for jets $p_{T}^{j}>30$ GeV and
$|\eta_{j}|<2.8$, and for the missing transverse momentum $p_{T}^{\rm
miss}>100$ GeV. It turns out the production cross sections with the precuts
are about $3.0$ fb and $3.8$ fb for BP 1 and 2 respectively, and dominated by
$Z$ and $Z^{\prime}$ mediated diagrams. We generate $10^{4}$ events for the
$pp\rightarrow W^{\prime p}W^{\prime m}j$ process and recast ATLAS mono-jet
search ATLAS:2020wzf  using MadAnalysis 5 Dumont:2014tja . The most sensitive
signal region is found to be in the window $E_{T}^{\rm miss}\in(700,800)$ GeV
(the signal region EM7 in Ref. ATLAS:2020wzf ). The 95% C.L. exclusion limits
on the production cross section are 400 fb and 680 fb for BP 1 and 2
respectively, which are much larger than the signal expected from the model.
Therefore, the LHC with luminosity of 139 fb-1 is not sensitive enough to
search for mono-jet events from this model. However, the model can be probed
by mono-jet searches at future hadron colliders such as the High-Luminosity
Large Hadron Collider (HL-LHC) Apollinari:2017cqg , the High-Energy Large
Hadron Collider (HE-LHC) Benedikt:2018ofy  and the Future Circular hadron-
hadron Collider (FCC-hh) Arkani-Hamed:2015vfh .

##  V Results

###  V.1  Methodology

The masses of the gauge bosons, $m_{Z}$, $m_{Z^{\prime}}$ and
$m_{A^{\prime}}$, as well as the constraints from Sec. III are calculated
through our own fortran codes, except for the Higgs invisible decay which can
be calculated together with the DM constraints. The DM constraints of Sec. IV,
in particular relic density, direct detection and indirect detection are
calculated using micrOMEGAs Belanger:2018ccd  and a set of model files
generated by FeynRules Alloul:2013bka . For the invisible decay branching
ratio of the Higgs, we take advantage of the use of CalcHEP Belyaev:2012qa
within micrOMEGAs to calculate the decay width along with the rest of the DM
constraints just mentioned.

All the points outside the theoretical constraints of Sec. III.1 are simply
rejected. By the same token, the dark photon constraints are used to reject
any parameter combination of $\varepsilon_{e,\mu}$-$m_{A^{\prime}}$ or
$\varepsilon_{e,\mu}$-$m_{Z^{\prime}}$ located inside the currently excluded
regions. The rest of the constraints in Sec. III are summed into a total
$\chi^{2}$ that also includes relic density and direct detection cross
section. In the case of direct detection experiments, where a limit is
reported at a 95% C.L. with null-signal assumption, we use a
$\chi_{\text{DD}}^{2}$ of the form

$\chi^{2}_{\text{DD}}=4.61\times\left(\frac{\sigma_{\text{theory}}}{\sigma_{\text{limit}}}\right)^{2},$
(73)

where the 4.61 factors allows $\chi_{\text{DD}}^{2}=4.61$ when we are exactly
at the 95% C.L. of this two-dimensional limit.333 For a one-tailed test, the
95% C.L. corresponds to $\Delta\chi^{2}=2.71$ and $\Delta\chi^{2}=4.61$ of a
Gaussian distribution in one and two dimensions, respectively. For a two-
tailed test, the same numbers correspond to the 90% C.L. limit. In mass ranges
where more than one limit exists we take the one with the largest
$\chi_{\text{DD}}^{2}$. Note that, due to ISV, the largest $\chi^{2}$ for
direct detection may not correspond to the experiment with the smallest cross
section. Since direct detection limits are reported assuming $f_{p}=f_{n}$ in
Eq. (71), it is possible for ISV ($f_{p}\neq f_{n}$) to produce some amount of
cancellation or enhancement of the limits depending on the atoms used in the
detector. Calculating the cross section in the way described in Sec. IV.2
allows us to account for ISV and the atoms used in different experiments.

In the case of Higgs invisible decay branching fraction, where a limit is
reported with a 95% C.L., the appropriate $\chi_{\text{inv}}^{2}$ is given by

$\chi^{2}_{\text{inv}}=2.71\times\left(\frac{{\rm BR}({h_{1}\rightarrow{\rm
inv}})}{0.13}\right)^{2},$ (74)

where, similarly to direct detection, the 2.71 factor allows for
$\chi^{2}_{\text{inv}}=2.71$ when our result is exactly at the reported 95%
C.L. in the one-dimensional case.

Parameter [units] | Range
---|---
$m_{h_{1}}$ [GeV] | [124.26 , 125.94]
$m_{h_{2}}$ [TeV] | [0.3 , 10]
$m_{D}$ [TeV] | [0.3 , 10]
$m_{H^{\pm}}$ [TeV] | [0.3 , 10]
$\theta_{1}$ [rad] | [$-\pi/2$ , $\pi/2$]
$\log_{10}(m_{W^{\prime}}/\text{GeV})$ | [$-3$ , 2]
$\log_{10}(M_{X}/\text{GeV})$ | [$-3$ , 2]
$\log_{10}(g_{H})$ | [$-6$ , 0]
$\log_{10}(g_{X})$ | [$-6$ , 0]
$m_{f^{\rm H}}$ [TeV] | 3 (fixed)
Table 3:  Ranges and values for the prior of the parameters used in this
analysis. All priors are taken as uniform inside their ranges for the
parameters listed in the table. While the $m_{h_{1}}$ range in this table
corresponds to the measured $\pm 6\sigma$ interval, emcee takes care of
sampling it according to the result of the total $\chi^{2}$.

To sample the parameter space we use the affine invariant Markov Chain Monte
Carlo (MCMC) ensemble sampler emcee ForemanMackey:2012ig  which presents
advantages such as fast calculation of parameter distributions in several
dimensions. The initial prior and ranges of each parameter are contained in
Table 3. In particular, the parameters $m_{W^{\prime}}$, $M_{X}$, $g_{H}$ and
$g_{X}$ are scanned in base-10 logarithmic scale. This is mostly because we
expect these parameters to be small but different from zero and that their
effects depend heavily on their orders of magnitude. For the masses of the
heavy fermions, we expect their contributions to be heavily suppressed by the
requirement that $m_{f^{\rm H}}=\mathcal{O}(\text{1~{}TeV})$. Therefore we
consider all of them degenerated with a mass of $m_{f^{\rm H}}=3$ TeV putting
them safely above any current search for heavy fermionic states. The rest of
the parameters are scanned with a uniform prior in linear scale. To guarantee
that our final distributions are independent of the initial points we perform
several small runs collecting $\mathcal{O}(10^{4})$ points for each run using
$\mathcal{O}(100)$ walkers. The initial points for the walkers are always
allowed by theoretical and dark photon constraints but otherwise random inside
the prior. After checking that the final distributions are consistent between
different runs, we perform a large scan with 300 walkers collecting 160,000
points after burn-in and thinning.

###  V.2  Numerical Results

We present the numerical results for visualization in Figs. 2–5. To follow the
discussion below more smoothly, we suggest our readers to read the captions of
these figures first and then view and compare them in parallel.

Figure 2:  The 1$\sigma$, 2$\sigma$ and 3$\sigma$ allowed regions (dark blue,
medium blue and light blue, respectively) projected on the planes
$(m_{W^{\prime}},g_{H})$ (left) and $(g_{H},g_{X})$ (right). The solid, dashed
and dotted black contours delimit the 1$\sigma$, 2$\sigma$ and 3$\sigma$
regions respectively. The labels on empty zones represent the most relevant
constraint in that region. The mono-jet benchmark points of Sec. IV.4 are
shown as a green down-triangle (BP 1) and a red up-triangle (BP 2).

The most notable feature of both panes of Fig. 2 is the band-shaped allowed
region. In the case of $(m_{W^{\prime}},g_{H})$ plane shown in the left pane,
the band is caused by the relation between relic density and cross section,
$\Omega h^{2}\propto 1/\langle\sigma v\rangle$. Considering that we have
$\sigma\propto g_{H}^{2}m_{W^{\prime}}^{2}/s^{2}$ from Eq. (IV.1), assuming
$s\sim 4\,m_{W^{\prime}}^{2}$ we have that
$g_{H}^{2}m_{W^{\prime}}^{2}/s^{2}\sim g_{H}^{2}/(16\,m_{W^{\prime}}^{2})$
resulting in $\Omega h^{2}\propto m_{W^{\prime}}^{2}/g_{H}^{2}$. This means
that to keep a constant relic density, $m_{W^{\prime}}$ and $g_{H}$ have to
keep a linear relationship as displayed in the left pane of Fig. 2. Deviations
from this band result in the relic density going either above or below the
value measured by the Planck satellite.

In the case of the right pane of Fig. 2, the band can be explained by the
possibility of having a resonant annihilation of $W^{\prime(p,m)}$ mediated by
the $Z^{\prime}$. First note that Eq. (23) directly relates $m_{A^{\prime}}$
and $M_{X}$ and, as can be seen in the right pane of Fig. 3, $m_{A^{\prime}}$
is required to be mostly below 0.2 GeV due to the LHCb results, thus limiting
also the size of $M_{X}$. Then, from Eq. (22) we know that the term with
$m_{W^{\prime}}^{2}$ factor dominates over the term with $M_{X}^{2}$. Finally,
resonant annihilation is achieved for $m_{Z^{\prime}}\approx
2\,m_{W^{\prime}}$ meaning $1+4g_{X}^{2}/g_{H}^{2}\approx 4$ or
$g_{X}^{2}/g_{H}^{2}\approx 3/4$, resulting in the band seen in the right pane
of Fig. 2. Again, large deviations from this band result in too much or not
enough annihilation to achieve the correct relic density. Due to this
$g_{X}$-$g_{H}$ correlation, the two-dimensional allowed regions in Fig. 4 for
all the distributions involving $g_{X}$ and $g_{H}$ have similar shapes. Here
it is important to mention that exact resonance,
$m_{Z^{\prime}}=2\,m_{W^{\prime}}$, would result in too much annihilation and,
therefore, in $\Omega h^{2}$ well below Planck’s measurement.

Besides their band-shaped tendency, both panes in Fig. 2 are bounded in their
top-right and bottom-left corners by the DM direct detection and dark photon
constraints, respectively. We know that the direct detection cross section
grows with $g_{H}^{2}$ as seen in Eq. (72), therefore, it is expected to see
it setting an upper bound on $g_{H}$. Furthermore, as can be seen in the left
pane of Fig. 3, direct detection experiments practically create a wall that
limits the size of $m_{W^{\prime}}$ from above. The effects of this limit are
reflected in the upper bound of $m_{W^{\prime}}$ in the left pane of Fig. 2.
In the case of the region disfavored by dark photon searches, this is mostly
due to the $\nu$-CAL I experiment limiting $\varepsilon$ from below as seen in
the right pane (olive green shaded zone) of Fig. 3. The $\varepsilon$ coupling
limit is passed to $g_{H}$ through the vectorial and axial couplings
$v_{f}^{A^{\prime}}$ and $a_{f}^{A^{\prime}}$ that depend on it, resulting on
the lower limit on $g_{H}$ that can be seen in both panes of Fig. 2.

Figure 3:  The 1$\sigma$ and 2$\sigma$ allowed contours projected on the DM
mass, $m_{W^{\prime}}$, vs. direct detection cross section (left) and dark
photon mass, $m_{A^{\prime}}$, vs. $\varepsilon$ coupling (right). The
experimental excluded regions used in this study are shown as solid colored
regions. Projected experimental limits are shown as dotted lines, with the
direction of the exclusion marked in gray. In the case of direct detection
(left), the limit set by the neutrino floor is shown as a dashed light orange
line.

In the left pane of Fig. 3 we show the allowed region projected on the
($m_{W^{\prime}}$, $\sigma_{W^{\prime}p}^{\rm SI}$) plane, where
$\sigma_{W^{\prime}p}^{\rm SI}$ is the cross section for spin-independent
scattering on a proton. The dark (light) blue shaded zone represents the
$1\sigma$ ($2\sigma$) allowed region. The current DM direct detection
measurements from CRESST III (green) Angloher:2017sxg , DarkSide-50 (orange)
Agnes:2018ves  and XENON1T (brown) Aprile:2019xxb  constrain the DM mass to
remain below $\sim 2$ GeV. A small part of the $2\sigma$ allowed region lies
below the neutrino floor (light orange), where the coherent neutrino-nucleus
scattering would dominate over any DM signal. Additionally, we show that
experiments in the near future such as NEWS-G Battaglieri:2017aum , SuperCDMS
Agnese:2016cpb  and CDEX Ma:2017nhc  can further probe our allowed parameter
space, in particular for $m_{W^{\prime}}\mathrel{\mathop{\kern 0.0pt\hbox
to0.0pt{\raise 0.86108pt\hbox{$>$}\hss}}\lower
3.87495pt\hbox{\kern-1.90002pt$\sim$}}0.3$ GeV with NEWS-G and down to
$\sigma_{W^{\prime}p}^{\rm SI}\sim 10^{-44}$ cm2 with SuperCDMS and CDEX.

The right pane in Fig. 3 shows the allowed region projected on the
($m_{A^{\prime}}$, $\varepsilon$) plane with the coupling
$\varepsilon\equiv\varepsilon_{\ell}$. Various experimental limits from dark
photon searches are displayed in color shaded zones including LHCb (green)
Aaij:2019bvg , BaBar (pink) Lees:2014xha , NA48 (purple) Batley:2015lha , NA64
(light brown) Banerjee:2018vgk , E141 (magenta) Riordan:1987aw  and $\nu$-CAL
I (light green) Blumlein:2011mv ; Blumlein:2013cua . The dilepton searches at
the LHCb, BaBar and NA48 put upper limits of $\varepsilon\lesssim 10^{-3}$ for
$m_{A^{\prime}}\mathrel{\mathop{\kern 0.0pt\hbox to0.0pt{\raise
0.86108pt\hbox{$>$}\hss}}\lower 3.87495pt\hbox{\kern-1.90002pt$\sim$}}0.03$
GeV, especially LHCb which sets a strong limit on $\varepsilon$ at $0.2$ GeV
$<m_{A^{\prime}}<0.5$ GeV causing a concave region in the $2\sigma$ allowed
region at this mass range. We note that this concave region due to LHCb
corresponds to the concave region at
$(m_{W}^{\prime},\sigma_{W^{\prime}p}^{\rm
SI})\sim(1\rm\,{GeV},10^{-42}\,\rm{cm^{2}})$ in the left pane of the same
figure. The LHCb long lived dark photon search constraints Aaij:2019bvg  are
also shown by the two isolated green shaded islands around $\epsilon$ equals
$2\times 10^{-5}$. On the other hand, the beam dump experiments NA64, E141 and
$\nu$-CAL I close the available space for smaller $\varepsilon$ and lighter
$m_{A^{\prime}}$ setting lower bounds of $m_{A^{\prime}}>0.02$ GeV and
$\varepsilon\mathrel{\mathop{\kern 0.0pt\hbox to0.0pt{\raise
0.86108pt\hbox{$>$}\hss}}\lower 3.87495pt\hbox{\kern-1.90002pt$\sim$}}2\times
10^{-5}$. The lower limit on $\varepsilon$ for $m_{A^{\prime}}>0.05$ GeV is
due to the DM relic density measured by the Planck experiment. Interestingly,
our final allowed region is located in the gap between the beam-dump and the
collider based experiments, an area of special interest for future dark photon
searches. For example, Belle-II Kou:2018nap  with a luminosity of $50\ {\rm
ab}^{-1}$ can probe $\varepsilon$ down to $2\times 10^{-4}$, the next upgrade
of NA64 NA64:2018  can cover $10^{-5}\mathrel{\mathop{\kern 0.0pt\hbox
to0.0pt{\raise 0.86108pt\hbox{$<$}\hss}}\lower
3.87495pt\hbox{\kern-1.90002pt$\sim$}}\varepsilon\mathrel{\mathop{\kern
0.0pt\hbox to0.0pt{\raise 0.86108pt\hbox{$<$}\hss}}\lower
3.87495pt\hbox{\kern-1.90002pt$\sim$}}10^{-3}$ and
$m_{A^{\prime}}\mathrel{\mathop{\kern 0.0pt\hbox to0.0pt{\raise
0.86108pt\hbox{$<$}\hss}}\lower 3.87495pt\hbox{\kern-1.90002pt$\sim$}}0.08$
GeV by reaching $\sim 5\times 10^{12}$ electrons-on-target (abbreviated by eot
in the figure) and Advanced WAKEfield Experiment (AWAKE) run 2
Caldwell:2018atq  can reach $m_{A^{\prime}}$ up to $0.15$ GeV with $10^{16}$
electrons-on-target with an energy of 50 GeV. These limits are shown
explicitly in the right pane of Fig. 3 as dotted lines with the side of the
exclusion in gray. In the future, with access to high energy electron-proton
colliders, AWAKE may reach 1 TeV for the electrons, extending $m_{A^{\prime}}$
up to 0.6 GeV Caldwell:2018atq  and dark photon searches at LHeC and FCC-he
DOnofrio:2019dcp  may even cover our entire allowed parameter space.

Figure 4:  Marginalized distributions in two dimensions for the parameters
scanned in this study. The 1$\sigma$, 2$\sigma$ and 3$\sigma$ regions are
marked in dark, light and lighter blue delimited by solid, dashed and dotted
lines, respectively.

We present two-dimensional projections of the allowed region for our BSM
parameters in Fig. 4. The dark, light and lighter blue zones indicate the
$1\sigma$, $2\sigma$ and $3\sigma$ allowed regions respectively. One important
thing to note is that the $m_{h_{2}}$, $m_{D}$ and $m_{H^{\pm}}$ masses show
an apparent upper bound at 1$\sigma$ and 2$\sigma$. This upper bound actually
depends on the maximum value chosen for the prior of these three parameters
and has no physical meaning. In the case of $m_{h_{2}}$, the apparent limit is
due to the reduced $\theta_{1}$ for large $m_{h_{2}}$ seen in the
$(m_{h_{2}},\theta_{1})$ subfigure while for $m_{D}$ and $m_{H^{\pm}}$ it is
due to their near degeneracy shown in the $(m_{D},m_{H^{\pm}})$ subfigure. We
have checked that changing the maximum scanned value for these three
parameters does not change the rest of the distributions, except for
$\theta_{1}$ where, understandably, larger $m_{h_{2}}$ sharpens the peak at
$\theta_{1}=0$ where $h_{1}=h_{\rm SM}$ exactly. The near degeneracy for
$m_{D}$ and $m_{H^{\pm}}$ is due to their comparably small mass squared
difference
$m_{D}^{2}-m_{H^{\pm}}^{2}=v^{2}(\lambda^{\prime}_{H\Phi}+\lambda^{\prime}_{H})/2$.
Given that $(\lambda^{\prime}_{H\Phi}+\lambda^{\prime}_{H})/2$ is bounded by
unitarity constraints, the mass squared difference is expected to remain
$\mathcal{O}(v^{2})$ or less, meaning that as $m_{D}$ and $m_{H^{\pm}}$ grow
away from $v$ their proportional difference rapidly grows smaller. Their near
degeneracy is also noticeable in all their two-dimensional distributions since
the distributions become nearly identical as $m_{D}$ and $m_{H^{\pm}}$ grow
larger.

Another interesting feature is that the charged Higgs mass distributions
reveal the presence of a lower limit around $m_{H^{\pm}}\approx 400$ GeV where
the contours show that the distribution falls rapidly. This is due to the
constraint from the Higgs decays into diphoton as shown in Eq. (52). Moreover,
due to the relation between $M_{X}$ and $m_{A^{\prime}}$ that can be inferred
from Eq. (23), the distributions in the $(M_{X},g_{X})$ and $(M_{X},g_{H})$
subfigures (second column from the right in (Fig. 4)) are close to the
distribution shown for the dark photon constraint in the right pane of Fig. 3.

Figure 5:  Marginalized distributions of the most relevant parameters used in
this study. Central values are indicated in the $x$-axis label. Since all the
distributions are normalized, the $y$-axis values are only for reference.

We present the marginalized one-dimensional distributions for the most
relevant parameters in Fig. 5 together with their central values and their
1$\sigma$ intervals. As mentioned before, there is a relation between $g_{X}$
and $g_{H}$ which is seen again in their marginalized distributions in the top
two panes since both peak close to $10^{-4}$ and have very similar 1$\sigma$
intervals. The distribution for $\log_{10}(M_{X}/\text{GeV})$ in the middle-
left pane shows the effects of the LHCb constraint around
$\log_{10}(M_{X}/\text{GeV})\approx 0$ and peaks just above $-1$. Because of
the precise measurements on the Higgs properties at the LHC, the mixing angle
between $h_{1}$ (identified to be the observed 125 GeV Higgs) and $h_{2}$,
$\theta_{1}$ is found to have a 1$\sigma$ interval of $(-9.8,9.5)\times
10^{-2}$ rad as shown in the middle-right pane of Fig. 5. For the distribution
of the DM mass $m_{W^{\prime}}$ in the bottom pane, we have the 1$\sigma$
interval for $\log_{10}(m_{W^{\prime}}/\text{GeV})$ between $-0.81$ and
$-0.04$ which corresponds to $m_{W^{\prime}}\in[0.15,0.91]$ GeV peaking at
$m_{W^{\prime}}\approx 0.39$ GeV.

After the above long discussions of our numerical analysis, perhaps a high
level summary is useful.

  * •

We have focused our numerical analysis on the parameter space of the model
that can lead to a sub-GeV $W^{\prime\,(p,m)}$ DM with a mass range of
MeV$-$GeV.

  * •

The viable domains of the parameter space are summarized in Fig. 4, where the
1–3$\sigma$ contours allowed by all existing experimental constraints are
shown for any combination of two of the 8 free parameters while the rest of
the parameters are marginalized. All the masses of the new heavy fermions in
the model, required by anomaly cancellation, are set at 3 TeV. For the two
cases of $(m_{W^{\prime}},g_{H})$ and $(g_{H},g_{X})$, more detailed
information on the boundaries of the contours due to the different constraints
imposed are exhibited in Fig. 2.

  * •

It is both clear and exciting to see from the allowed $1\sigma$ and $2\sigma$
contours on the plane of $m_{W^{\prime}}$ versus $\sigma^{\rm
SI}_{W^{\prime}p}$ displayed in the left pane of Fig. 3 that future
experiments like NEWS-G and SuperCDMS can put more stringent constraints on
the model. Indeed, about 1/3 of the current viable parameter space in the
$(m_{W^{\prime}},\sigma^{\rm SI}_{W^{\prime}p})$ plane would be facing
challenge. On the other hand, the region where the projected CDEX sensitivity
can reach is already disfavored in G2HDM. Note that a small portion of the
$2\sigma$ allowed parameter space in the $(m_{W^{\prime}},\sigma^{\rm
SI}_{W^{\prime}p})$ plane is overlapping with the neutrino floor.

  * •

For the allowed contours of dark photon coupling $\varepsilon$ and mass
$m_{A^{\prime}}$ exhibited in the right pane of Fig. 3, future upgraded NA64
experiment with $5\times 10^{12}$ electrons-on-target, proof-of-principle
experiment AWAKE run 2 with $10^{16}$ electrons-on-target, and next generation
B-factory experiment Belle II would probe more than 1/2 of the current viable
domain in the $(m_{A^{\prime}},\varepsilon)$ plane.

  * •

The correlation between the DM and dark photon physics and their constraints
in G2HDM exhibited in the left and right panes of Fig. 3 is rather novel and
interesting. Indeed any non-abelian vector DM is likely accompanied by at
least one extra neutral gauge boson which can play the role of $Z^{\prime}$ or
$A^{\prime}$. Thus considerations of experimental constraints for both DM and
dark photon physics must be taken into account mandatory. Although our
analysis is performed in the context of a specific model, we expect some of
the features obtained in this work may be generic for any low mass non-abelian
vector DM with a dominated vector portal of neutral gauge boson communicating
to the SM sector. Perhaps we may be entertained by nature revealing to us not
only a low mass dark matter $W^{\prime}$ but also a dark photon and a dark
$Z^{\prime}$ with nearby masses. These are the three gauge bosons associated
with the dark $SU(2)_{H}\times U(1)_{X}$ sector, in mirror with the visible SM
gauge group $SU(2)_{L}\times U(1)_{Y}$. Due to the minuscule couplings of
$g_{H}$ and $g_{X}$, direct detections of the dark matter, dark photon and
dark $Z^{\prime}$ (as well as other new particles introduced in G2HDM) at
colliders would belong to the lifetime frontier and high luminosity/energy
frontier at the upgrade of LHC and future colliders.

##  VI Conclusion

In this work, we studied a simplified version of the G2HDM Huang:2015wts . In
previous works, it was demonstrated that the original G2HDM successfully
explains dark matter while keeping a parameter space consistent with
theoretical and experimental expectations Arhrib:2018sbz ; Huang:2019obt ;
Chen:2019pnt . The simplifications considered in this study reduced the size
of the scalar potential and the parameter space by removing the scalar
$SU(2)_{H}$ triplet. Thankfully, the absence of this triplet does not affect
the $h$-parity and another dark sector particle, the $W^{\prime(p,m)}$, steps
up as a DM candidate. For the properties of $W^{\prime(p,m)}$ we settled on
exploring the sub-GeV mass range. Interestingly, in this region the extra
vector states $Z^{\prime}$ and $A^{\prime}$ happen to have relevant roles as
mediators for the annihilation and direct detection processes of the DM. In
particular, reaching the appropriate amount of annihilation to have the
correct relic density would have been difficult if not impossible without the
channel mediated by the $Z^{\prime}$.

We started with the usual theoretical checks on the scalar potential ensuring
that the minimum is stable and that the couplings remain unitary at tree
level. Given that the LHC has been closing in on the detailed properties of
the 125 GeV Higgs, we check that our scalar sector provides a particle,
$h_{1}$, that matches the mass and decays that have been measured. The same
can be said in the case of the gauge sector and the $Z$ with its properties
already very well measured at LEP. In the case of the light $Z^{\prime}$ and
$A^{\prime}$ we constrained their interactions with SM leptons by checking
against the regions excluded by dark photon searches in LHCb, BaBar, NA48,
NA64, E141 and $\nu$-CAL I. Finally, we required our DM candidate,
$W^{\prime(p,m)}$, to have the correct relic density measured by Planck
satellite with a spin-independent direct detection cross section below the
limits found by CRESST-III, DarkSide-50 and XENON1T.

In our numerical analysis, we uncovered some interesting features of the
parameter space, such as a correlation between the couplings $g_{H}$ and
$g_{X}$ with superweak size $O(10^{-5}-10^{-3})$ and that most of our results
lie inside the gap between various dark photon explorations. The latter
becomes more important when we consider that this gap is projected to be
further explored in the future with the upgrades to NA64 and AWAKE and the
B-factory Belle II. This would reduce the allowed parameter space nearly by
half, and therefore remove a good portion of our allowed parameter space in
the lighter $A^{\prime}$ region. Moreover, future DM direct detection
experiments like NEWS-G and SuperCDMS may reduce the parameter space by
exploring the regions with heavier $W^{\prime(p,m)}$ and larger cross section.

To summarize, we found that the simplified G2HDM developed in this work
provides a viable vector DM candidate with mass down to $\mathcal{O}(10^{-2})$
GeV. All the predictions in the scalar and gauge sectors are in good agreement
with current observations. Importantly, both new vector states, $Z^{\prime}$
and $A^{\prime}$, play key roles for DM observables. Besides the possibility
of detecting the $W^{\prime(p,m)}$ in DM direct detection experiments, the
dark photon, $A^{\prime}$, is predicted to be well positioned for future
observations that may reach $m_{A^{\prime}}\sim 0.1$ GeV. This work
demonstrates that the G2HDM is not only a successful and competitive dark
matter formulation but can also serve as a starting point with diverse
exploration possibilities.

## Acknowledgments

The analysis presented here was done using the resources of the high-
performance T3 Cluster at the Institute of Physics, Academia Sinica. This work
is supported in part by the Ministry of Science and Technology of Taiwan under
Grants No. 108-2112-M-001-046 (TCY) and No. 109-2811-M-001-595 (RR) and by
National Natural Science Foundation of China under Grants No. 11775109 (VQT).
VQT would like to thank the Institute of Physics, Academia Sinica, Taiwan for
its hospitality during this work.

##  Appendix A Mass Spectra of Goldstone Bosons and Gauge Fixings in General
Renormalizable Gauge

In Sec. II.3 as well as in all previous works, we have derived the scalar and
vector boson mass spectra in the ’t Hooft-Landau gauge in which all the gauge
parameter parameters $\xi$s were set to zero. This hides away the issue of
gauge dependence and all the Goldstone bosons are massless. Since the physical
dark Higgs $D$ is a linear combination of $H_{2}^{0*}$ in the inert doublet
$H_{2}$ and Goldstone boson $G^{p}_{H}$ in the hidden doublet $\Phi_{H}$, one
might wonder what would happen in the general renormalizable gauge. We take
the opportunity here to examine this question more careful and discuss the
related issue of gauge mixings in G2HDM.

###  A.1 $G^{\pm}$ and $G^{(p,m)}_{H}$

First $G^{\pm}$. This is the same as in SM, which after SSB the covariant
kinetic term of the Higgs field $H_{1}$ contains the following mixing term

$-\frac{i}{2}gvW^{+}_{\mu}\partial^{\mu}G^{-}+{\rm H.c.}$ (75)

In the general renormalizable gauge, one introduces the following gauge fixing
term

$-\frac{1}{\xi_{W}}\Bigl{|}\partial^{\mu}W_{\mu}^{+}-\frac{i}{2}\xi_{W}gvG^{+}\Bigr{|}^{2}\;,$
(76)

where $\xi_{W}$ is an arbitrary gauge parameter. Expanding out Eq. (76), up to
a total derivative, it will cancel the mixing term in Eq. (75). In addition,
it will modify the $W^{\pm}$-boson propagator to be $\xi_{W}$-dependent whose
form is well known in the literature, and the Goldstone boson $G^{\pm}$
(absorbed by the longitudinal component of $W^{\pm}$) will develop a mass
equals to $gv\sqrt{\xi_{W}}/2$.

The story of $G^{(p,m)}_{H}$ in G2HDM is a little bit more interesting. The
mixing term contains two contributions coming from the covariant kinetic terms
of the $H_{2}$ and $\Phi_{H}$ fields

$-\frac{i}{2}g_{H}W^{p}_{\mu}\left(v_{\Phi}\partial^{\mu}G^{m}_{H}-v\partial^{\mu}H^{0}_{2}\right)\;,$
(77)

which indicates that the physical Goldstone field $\tilde{G}^{p}_{H}$
($\tilde{G}^{m}_{H}$) is actually a linear combination of $H^{0*}_{2}$
($H^{0}_{2}$ ) and $G^{p}_{H}$ ($G^{m}_{H}$), viz.

$\tilde{G}^{p}_{H}=\frac{1}{\sqrt{v^{2}+v_{\Phi}^{2}}}\left(v_{\Phi}G^{p}_{H}-vH^{0*}_{2}\right)\;,\;\;\tilde{G}^{m}_{H}=\frac{1}{\sqrt{v^{2}+v_{\Phi}^{2}}}\left(v_{\Phi}G^{m}_{H}-vH^{0}_{2}\right)\;.$
(78)

${\tilde{G}}^{(p,m)}_{H}$ are absorbed by the longitudinal components of
$W^{\prime(p,m)}$. The other physical orthogonal combination is the complex
dark Higgs $D$

$D=\frac{1}{\sqrt{v^{2}+v_{\Phi}^{2}}}\left(vG^{p}_{H}+v_{\Phi}H^{0*}_{2}\right)\;,\;\;D^{*}=\frac{1}{\sqrt{v^{2}+v_{\Phi}^{2}}}\left(vG^{m}_{H}+v_{\Phi}H^{0}_{2}\right)\;.$
(79)

In analogous with Eqs. (75) and (76), to cancel the mixing term in Eq. (77)
one introduces the following general gauge fixing term

$-\frac{1}{\xi_{W^{\prime}}}\Bigl{|}\partial^{\mu}W_{\mu}^{p}-\frac{i}{2}\xi_{W^{\prime}}g_{H}\tilde{v}\tilde{G}^{p}_{H}\Bigr{|}^{2}\;,$
(80)

where $\xi_{W^{\prime}}$ is an arbitrary gauge parameter and
$\tilde{v}=\sqrt{v^{2}+v^{2}_{\Phi}}$. As in the case of $W^{\pm}$, the
propagator of $W^{\prime(p,m)}$ is then get modified under this general
renormalizable gauge. The mass matrix ${\mathcal{M}}^{2}_{S^{\prime}}$ in Eq.
(12) is also modified as

${\mathcal{M}}_{S^{\prime}}^{2}=\begin{pmatrix}\frac{1}{2}\lambda^{\prime}_{H\Phi}v^{2}+\frac{1}{4}\xi_{W^{\prime}}g_{H}^{2}v_{\Phi}^{2}&\frac{1}{2}\lambda^{\prime}_{H\Phi}vv_{\Phi}-\frac{1}{4}\xi_{W^{\prime}}g_{H}^{2}vv_{\Phi}\\\
\frac{1}{2}\lambda^{\prime}_{H\Phi}vv_{\Phi}-\frac{1}{4}\xi_{W^{\prime}}g_{H}^{2}vv_{\Phi}&\frac{1}{2}\lambda^{\prime}_{H\Phi}v_{\Phi}^{2}+\frac{1}{4}\xi_{W^{\prime}}g_{H}^{2}v^{2}\end{pmatrix}\;,$
(81)

which has two eigenvalues $g_{H}^{2}\xi_{W^{\prime}}\tilde{v}^{2}/4$ and
$\lambda^{\prime}_{H\Phi}\tilde{v}^{2}/2$. The first one is the mass-squared
of the Goldstone boson $\tilde{G}^{(p,m)}_{H}$ in the general gauge, while the
second one is the mass-squared of the complex dark Higgs $D$ which is the same
as Eq. (16) previously derived from the ’t Hooft-Landau gauge as it should!

###  A.2 $G^{0}$, $G^{0}_{H}$ and $\mathcal{S}$

Similarly for the neutral gauge bosons, after SSB, we have the following
mixing terms from the covariant kinetic terms of the $H_{1}$ and $\Phi_{H}$
fields

$\displaystyle-$
$\displaystyle\\!\\!\frac{1}{2}v\left(g^{\prime}B^{\mu}-gW^{3\mu}+g_{H}W^{\prime
3\mu}+2g_{X}X^{\mu}\right)\partial_{\mu}G^{0}$ (82) $\displaystyle+$
$\displaystyle\\!\\!\frac{1}{2}v_{\Phi}\left(g_{H}W^{\prime
3\mu}-2g_{X}X^{\mu}\right)\partial_{\mu}G^{0}_{H}\;.$

For the $U(1)_{X}$ gauge field, we also use the Stueckelberg mechanism to
provide a mass. The Stueckelberg Lagrangian is

$\displaystyle{\mathcal{L}}_{\rm St}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(\partial^{\mu}\mathcal{S}+M_{X}X^{\mu}\right)^{2}\;,$
(83) $\displaystyle=$
$\displaystyle\frac{1}{2}\left(\partial^{\mu}\mathcal{S}\right)^{2}+\frac{1}{2}M_{X}^{2}X^{\mu}X_{\mu}+M_{X}X^{\mu}\partial_{\mu}\mathcal{S}\;.$

Here $\mathcal{S}$ and $M_{X}$ are the Stueckelberg field and mass
respectively. The middle term gives $X$ a mass $M_{X}$, while the last term
indicates the mixing of $\mathcal{S}$ with the longitudinal component of
$X^{\mu}$. Note that $\mathcal{S}$ is massless.

Recall that the first two terms in Eq. (82) are SM-like and they can be
combined as

$(g^{\prime}B_{\mu}-gW^{3}_{\mu})=-(g^{\prime}\sin\theta_{W}+g\cos\theta_{W})Z^{\rm
SM}_{\mu}\;.$ (84)

Thus the total mixing terms including both Eq. (82) and the last term in Eq.
(83) are

$\displaystyle+$
$\displaystyle\\!\\!\frac{1}{2}v\left[(g^{\prime}\sin\theta_{W}+g\cos\theta_{W})Z^{{\rm
SM}\,\mu}-g_{H}W^{\prime 3\mu}-2g_{X}X^{\mu}\right]\partial_{\mu}G^{0}$ (85)
$\displaystyle+$ $\displaystyle\\!\\!\frac{1}{2}v_{\Phi}\left(g_{H}W^{\prime
3\mu}-2g_{X}X^{\mu}\right)\partial_{\mu}G^{0}_{H}+M_{X}X^{\mu}\partial_{\mu}\mathcal{S}\;.$

Note that, as one would expect, the photon field $A^{\mu}$ doesn’t enter in
Eq. (85) which is coming entirely from SSB and Stueckelberg mechanism. The
photon field remains massless, it has no associated Goldstone boson and its
general gauge fixing term is simply given by
$-\left(\partial_{\mu}A^{\mu}\right)^{2}/2\xi_{\gamma}$ as in the SM case. So
does the massless gluon field, whose gauge fixing term is
$-\left(\partial_{\mu}A^{a\mu}\right)^{2}/2\xi_{g}$ where $a=1,\ldots,8$ is
the adjoint index of the color group $SU(3)_{C}$.

As mentioned in the text, the neutral vector gauge bosons $Z^{\rm SM}$,
$W^{\prime 3}$ and $X$ are in general mixed together according to

$\begin{pmatrix}Z^{\rm SM}\\\ W^{\prime 3}\\\ X\end{pmatrix}={\cal
O}\cdot\begin{pmatrix}Z_{1}\\\ Z_{2}\\\ Z_{3}\end{pmatrix}\;,$ (86)

where $\cal O$ is an orthogonal matrix. In our numerical scan of the parameter
space in this work, we have $Z_{1}=Z\simeq Z^{\rm SM}$ which is very close to
the SM $Z$-boson, and $Z_{2}=Z^{\prime}$ and $Z_{3}=A^{\prime}$, both of which
are lighter than the $Z$. We will use $Z_{i}$ in this appendix instead of $Z$,
$Z^{\prime}$ and $A^{\prime}$.

In terms of the physical fields $Z_{i}$, Eq. (85) becomes

$\displaystyle+$
$\displaystyle\\!\\!\frac{1}{2}v\left[(g^{\prime}\sin\theta_{W}+g\cos\theta_{W}){\mathcal{O}}_{1i}Z_{i}^{\mu}-g_{H}{\mathcal{O}}_{2i}Z_{i}^{\mu}-2g_{X}{\mathcal{O}}_{3i}Z_{i}^{\mu}\right]\partial_{\mu}G^{0}$
(87) $\displaystyle+$
$\displaystyle\\!\\!\frac{1}{2}v_{\Phi}\left(g_{H}{\mathcal{O}}_{2i}Z_{i}^{\mu}-2g_{X}{\mathcal{O}}_{3i}Z_{i}^{\mu}\right)\partial_{\mu}G^{0}_{H}+M_{X}{\mathcal{O}}_{3i}Z_{i}^{\mu}\partial_{\mu}\mathcal{S}$
$\displaystyle\equiv$
$\displaystyle\sum_{i=1}^{3}Z_{i}^{\mu}\partial_{\mu}G_{i}\;.$

Here we have defined three physical Goldstone fields absorbed by the
longitudinal components of the three physical $Z_{i}$ as

$G_{i}={\widetilde{C}}_{i1}G^{0}+{\widetilde{C}}_{i2}G^{0}_{H}+{\widetilde{C}}_{i3}\mathcal{S}\;\;\;\;(i=1,2,3)\;,$
(88)

where the coefficients ${\widetilde{C}}_{ij}$ can be read off from the first
two lines in Eq. (87), namely

$\displaystyle{\widetilde{C}}_{i1}$ $\displaystyle=$
$\displaystyle\frac{1}{2}v\left[(g^{\prime}\sin\theta_{W}+g\cos\theta_{W}){\mathcal{O}}_{1i}-g_{H}{\mathcal{O}}_{2i}-2g_{X}{\mathcal{O}}_{3i}\right]\;,$
(89) $\displaystyle{\widetilde{C}}_{i2}$ $\displaystyle=$
$\displaystyle\frac{1}{2}v_{\Phi}\left(g_{H}{\mathcal{O}}_{2i}-g_{X}{\mathcal{O}}_{3i}\right)\;,$
(90) $\displaystyle{\widetilde{C}}_{i3}$ $\displaystyle=$ $\displaystyle
M_{X}{\mathcal{O}}_{3i}\;.$ (91)

To cancel the mixing term in the last expression in Eq. (87), we need the
following gauge fixing term,

$-\frac{1}{2\xi_{i}}\left(\partial_{\mu}Z^{\mu}_{i}-\xi_{i}G_{i}\right)^{2}\;,$
(92)

where $\xi_{i}$ are three arbitrary gauge parameters. This gauge fixing term
will not only modify the propagator of $Z_{i}$ but also induce a $3\times 3$
mass mixing matrix ${\mathcal{M}}^{2}_{G}$ among the three fields
$\\{G^{0},G_{H}^{0},\mathcal{S}\\}$ to be $\xi_{i}$-dependence, namely

${\mathcal{M}}^{2}_{G}=\begin{pmatrix}{\mathcal{M}}^{2}_{11}&{\mathcal{M}}^{2}_{12}&{\mathcal{M}}^{2}_{13}\\\
{\mathcal{M}}^{2}_{21}&{\mathcal{M}}^{2}_{22}&{\mathcal{M}}^{2}_{23}\\\
{\mathcal{M}}^{2}_{31}&{\mathcal{M}}^{2}_{32}&{\mathcal{M}}^{2}_{33}\end{pmatrix}\;\;\;,$
(93)

with the following matrix elements

$\begin{matrix}{\mathcal{M}}^{2}_{11}=\xi_{1}{\widetilde{C}}_{11}^{2}+\xi_{2}{\widetilde{C}}_{21}^{2}+\xi_{3}{\widetilde{C}}_{31}^{2}\;,\\\
{\mathcal{M}}^{2}_{12}=\xi_{1}{\widetilde{C}}_{11}{\widetilde{C}}_{12}+\xi_{2}{\widetilde{C}}_{21}{\widetilde{C}}_{22}+\xi_{3}{\widetilde{C}}_{31}{\widetilde{C}}_{32}={\mathcal{M}}^{2}_{21}\;,\\\
{\mathcal{M}}^{2}_{13}=\xi_{1}{\widetilde{C}}_{11}{\widetilde{C}}_{13}+\xi_{2}{\widetilde{C}}_{21}{\widetilde{C}}_{23}+\xi_{3}{\widetilde{C}}_{31}{\widetilde{C}}_{33}={\mathcal{M}}^{2}_{31}\;,\\\
{\mathcal{M}}^{2}_{22}=\xi_{1}{\widetilde{C}}_{12}^{2}+\xi_{2}{\widetilde{C}}_{22}^{2}+\xi_{3}{\widetilde{C}}_{32}^{2}\;,\\\
{\mathcal{M}}^{2}_{23}=\xi_{1}{\widetilde{C}}_{12}{\widetilde{C}}_{13}+\xi_{2}{\widetilde{C}}_{22}{\widetilde{C}}_{23}+\xi_{3}{\widetilde{C}}_{32}{\widetilde{C}}_{33}={\mathcal{M}}^{2}_{32}\;,\\\
{\mathcal{M}}^{2}_{33}=\xi_{1}{\widetilde{C}}_{13}^{2}+\xi_{2}{\widetilde{C}}_{23}^{2}+\xi_{3}{\widetilde{C}}_{33}^{2}\;.\end{matrix}$
(94)

While the eigenvalues of ${\mathcal{M}}^{2}_{G}$ are complicated functions of
the gauge fixing parameters, its determinant takes a simple form,

$\displaystyle{\rm Det}{\mathcal{M}}^{2}_{G}$ $\displaystyle=$
$\displaystyle\xi_{1}\xi_{2}\xi_{3}\biggl{(}{\widetilde{C}}_{11}{\widetilde{C}}_{23}{\widetilde{C}}_{32}+{\widetilde{C}}_{12}{\widetilde{C}}_{21}{\widetilde{C}}_{33}+{\widetilde{C}}_{13}{\widetilde{C}}_{22}{\widetilde{C}}_{31}\biggr{.}$
(95)
$\displaystyle\;\;\;\;\;\;\;\;\,\biggl{.}-{\widetilde{C}}_{13}{\widetilde{C}}_{21}{\widetilde{C}}_{32}-{\widetilde{C}}_{11}{\widetilde{C}}_{22}{\widetilde{C}}_{33}-{\widetilde{C}}_{12}{\widetilde{C}}_{23}{\widetilde{C}}_{31}\biggr{)}^{2}\;.$

Thus as any one of the gauge fixing parameters $\xi_{i}\rightarrow 0$, the
corresponding Goldstone boson $G_{i}$ has vanishing mass, reproducing the
result of ’t Hooft-Landau gauge.

This completes our discussions of all the gauge fixing terms in G2HDM. With
these gauge fixing terms at hand, one can straightforwardly obtain the gauge
fixing functions and the procedure of Faddeev-Popov quantization can be
proceeded as usual. Of course physical observables if computed correctly
should be independent of all the otherwise arbitrary gauge fixing parameters
$\xi$s discussed in this Appendix! In an ideal world, one would compute things
using arbitrary $\xi$s and show the dependence of $\xi$s is completely dropped
out at the end for any physical observable. In practice, things are hardly get
done that way.

## References

  * (1) V. Silveira and A. Zee, “Scalar Phantoms,” Phys. Lett. B 161, 136-140 (1985)
  * (2) N. G. Deshpande and E. Ma, “Pattern of Symmetry Breaking with Two Higgs Doublets,” Phys. Rev. D 18, 2574 (1978)
  * (3) E. Ma, “Verifiable radiative seesaw mechanism of neutrino mass and dark matter,” Phys. Rev. D 73, 077301 (2006) [arXiv:hep-ph/0601225 [hep-ph]].
  * (4) A. H. Chamseddine, R. L. Arnowitt and P. Nath, “Locally Supersymmetric Grand Unification,” Phys. Rev. Lett. 49, 970 (1982)
  * (5) P. Nath, R. L. Arnowitt and A. H. Chamseddine, “Gravity Induced Symmetry Breaking and Ground State of Local Supersymmetric {GUTs},” Phys. Lett. B 121, 33-36 (1983)
  * (6) H. P. Nilles, “Supersymmetry, Supergravity and Particle Physics,” Phys. Rept. 110, 1-162 (1984)
  * (7) H. C. Cheng and I. Low, “TeV symmetry and the little hierarchy problem,” JHEP 09, 051 (2003) [arXiv:hep-ph/0308199 [hep-ph]].
  * (8) H. C. Cheng and I. Low, “Little hierarchy, little Higgses, and a little symmetry,” JHEP 08, 061 (2004) [arXiv:hep-ph/0405243 [hep-ph]].
  * (9) I. Low, “T parity and the littlest Higgs,” JHEP 10, 067 (2004) [arXiv:hep-ph/0409025 [hep-ph]].
  * (10) J. McDonald, “Gauge singlet scalars as cold dark matter,” Phys. Rev. D 50, 3637-3649 (1994) [arXiv:hep-ph/0702143 [hep-ph]].
  * (11) C. P. Burgess, M. Pospelov and T. ter Veldhuis, “The Minimal model of nonbaryonic dark matter: A Singlet scalar,” Nucl. Phys. B 619, 709-728 (2001) [arXiv:hep-ph/0011335 [hep-ph]].
  * (12) 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]].
  * (13) R. Barbieri, L. J. Hall and V. S. Rychkov, “Improved naturalness with a heavy Higgs: An Alternative road to LHC physics,” Phys. Rev. D 74, 015007 (2006) [arXiv:hep-ph/0603188 [hep-ph]].
  * (14) L. Lopez Honorez, E. Nezri, J. F. Oliver and M. H. G. Tytgat, “The Inert Doublet Model: An Archetype for Dark Matter,” JCAP 02, 028 (2007) [arXiv:hep-ph/0612275 [hep-ph]].
  * (15) 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]].
  * (16) A. Belyaev, G. Cacciapaglia, I. P. Ivanov, F. Rojas-Abatte and M. Thomas, “Anatomy of the Inert Two Higgs Doublet Model in the light of the LHC and non-LHC Dark Matter Searches,” Phys. Rev. D 97, no.3, 035011 (2018) [arXiv:1612.00511 [hep-ph]].
  * (17) Y. L. S. Tsai, V. Tran and C. T. Lu, “Confronting dark matter co-annihilation of Inert two Higgs Doublet Model with a compressed mass spectrum,” JHEP 06, 033 (2020) [arXiv:1912.08875 [hep-ph]].
  * (18) S. Fabian, F. Goertz and Y. Jiang, “Dark Matter and Nature of Electroweak Phase Transition with an Inert Doublet,” [arXiv:2012.12847 [hep-ph]].
  * (19) S. Akula, B. Altunkaynak, D. Feldman, P. Nath and G. Peim, “Higgs Boson Mass Predictions in SUGRA Unification, Recent LHC-7 Results, and Dark Matter,” Phys. Rev. D 85, 075001 (2012) [arXiv:1112.3645 [hep-ph]].
  * (20) R. L. Arnowitt and P. Nath, “Predictions of neutralino dark matter event rates in minimal supergravity unification,” Phys. Rev. D 54, 2374-2384 (1996) [arXiv:hep-ph/9509260 [hep-ph]].
  * (21) G. Jungman, M. Kamionkowski and K. Griest, “Supersymmetric dark matter,” Phys. Rept. 267, 195-373 (1996) [arXiv:hep-ph/9506380 [hep-ph]].
  * (22) C. S. Chen, K. Cheung and T. C. Yuan, “Novel Collider Signature for Little Higgs Dark Matter Models,” Phys. Lett. B 644, 158-164 (2007) [arXiv:hep-ph/0605314 [hep-ph]].
  * (23) C. D. Carone and R. Ramos, “Classical scale-invariance, the electroweak scale and vector dark matter,” Phys. Rev. D 88, 055020 (2013) [arXiv:1307.8428 [hep-ph]].
  * (24) H. Davoudiasl and I. M. Lewis, “Dark Matter from Hidden Forces,” Phys. Rev. D 89, no.5, 055026 (2014) [arXiv:1309.6640 [hep-ph]].
  * (25) B. Barman, S. Bhattacharya, S. K. Patra and J. Chakrabortty, “Non-Abelian Vector Boson Dark Matter, its Unified Route and signatures at the LHC,” JCAP 12, 021 (2017) [arXiv:1704.04945 [hep-ph]].
  * (26) B. Barman, S. Bhattacharya and M. Zakeri, “Non-Abelian Vector Boson as FIMP Dark Matter,” JCAP 02, 029 (2020) [arXiv:1905.07236 [hep-ph]].
  * (27) T. Abe, M. Fujiwara, J. Hisano and K. Matsushita, “A model of electroweakly interacting non-abelian vector dark matter,” JHEP 07, 136 (2020) [arXiv:2004.00884 [hep-ph]].
  * (28) Z. Hu, C. Cai, Y. L. Tang, Z. H. Yu and H. H. Zhang, JHEP 07, 089 (2021) [arXiv:2103.00220 [hep-ph]].
  * (29) Q. H. Cao, C. R. Chen, C. S. Li and H. Zhang, “Effective Dark Matter Model: Relic density, CDMS II, Fermi LAT and LHC,” JHEP 08, 018 (2011) [arXiv:0912.4511 [hep-ph]].
  * (30) J. Goodman, M. Ibe, A. Rajaraman, W. Shepherd, T. M. P. Tait and H. B. Yu, “Constraints on Dark Matter from Colliders,” Phys. Rev. D 82, 116010 (2010) [arXiv:1008.1783 [hep-ph]].
  * (31) K. Cheung, P. Y. Tseng, Y. L. S. Tsai and T. C. Yuan, “Global Constraints on Effective Dark Matter Interactions: Relic Density, Direct Detection, Indirect Detection, and Collider,” JCAP 05, 001 (2012) [arXiv:1201.3402 [hep-ph]].
  * (32) J. Goodman, M. Ibe, A. Rajaraman, W. Shepherd, T. M. P. Tait and H. B. Yu, “Constraints on Light Majorana dark Matter from Colliders,” Phys. Lett. B 695, 185-188 (2011) [arXiv:1005.1286 [hep-ph]].
  * (33) J. Goodman, M. Ibe, A. Rajaraman, W. Shepherd, T. M. P. Tait and H. B. Yu, “Gamma Ray Line Constraints on Effective Theories of Dark Matter,” Nucl. Phys. B 844, 55-68 (2011) [arXiv:1009.0008 [hep-ph]].
  * (34) K. Cheung, P. Y. Tseng and T. C. Yuan, “Cosmic Antiproton Constraints on Effective Interactions of the Dark Matter,” JCAP 01, 004 (2011) [arXiv:1011.2310 [hep-ph]].
  * (35) K. Cheung, P. Y. Tseng and T. C. Yuan, “Gamma-ray Constraints on Effective Interactions of the Dark Matter,” JCAP 06, 023 (2011) [arXiv:1104.5329 [hep-ph]].
  * (36) W. C. Huang, K. W. Ng and T. C. Yuan, “Circularly Polarized Gamma Rays in Effective Dark Matter Theory,” Phys. Lett. B 800, 135104 (2020) [arXiv:1907.02402 [hep-ph]].
  * (37) T. Nath Maity and F. S. Queiroz, [arXiv:2104.02700 [hep-ph]].
  * (38) 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]].
  * (39) C. R. Chen, Y. X. Lin, C. S. Nugroho, R. Ramos, Y. L. S. Tsai and T. C. Yuan, “Complex scalar dark matter in the gauged two-Higgs-doublet model,” Phys. Rev. D 101, no.3, 035037 (2020) [arXiv:1910.13138 [hep-ph]].
  * (40) R. N. Mohapatra and G. Senjanovic, “Neutrino Mass and Spontaneous Parity Nonconservation,” Phys. Rev. Lett. 44, 912 (1980)
  * (41) R. N. Mohapatra and G. Senjanovic, “Neutrino Masses and Mixings in Gauge Models with Spontaneous Parity Violation,” Phys. Rev. D 23, 165 (1981)
  * (42) For a general review of 2HDM, see for example 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]].
  * (43) A. Arhrib, W. C. Huang, R. Ramos, Y. L. S. Tsai and T. C. Yuan, “Consistency of a gauged two-Higgs-doublet model: Scalar sector,” Phys. Rev. D 98, no.9, 095006 (2018) [arXiv:1806.05632 [hep-ph]].
  * (44) C. T. Huang, R. Ramos, V. Q. Tran, Y. L. S. Tsai and T. C. Yuan, “Consistency of Gauged Two Higgs Doublet Model: Gauge Sector,” JHEP 1909 (2019) 048 [arXiv:1905.02396 [hep-ph]].
  * (45) G. ’t Hooft, “Magnetic Monopoles in Unified Gauge Theories,” Nucl. Phys. B 79, 276-284 (1974)
  * (46) A. M. Polyakov, “Particle Spectrum in the Quantum Field Theory,” JETP Lett. 20, 194-195 (1974) PRINT-74-1566 (LANDAU-INST).
  * (47) B. Julia and A. Zee, “Poles with Both Magnetic and Electric Charges in Nonabelian Gauge Theory,” Phys. Rev. D 11, 2227-2232 (1975)
  * (48) S. Baek, P. Ko and W. I. Park, “Hidden sector monopole, vector dark matter and dark radiation with Higgs portal,” JCAP 10, 067 (2014) [arXiv:1311.1035 [hep-ph]].
  * (49) N. Aghanim et al. [Planck], “Planck 2018 results. VI. Cosmological parameters,” Astron. Astrophys. 641, A6 (2020) [arXiv:1807.06209 [astro-ph.CO]].
  * (50) G. Angloher et al. [CRESST], “Results on MeV-scale dark matter from a gram-scale cryogenic calorimeter operated above ground,” Eur. Phys. J. C 77, no.9, 637 (2017) [arXiv:1707.06749 [astro-ph.CO]].
  * (51) P. Agnes et al. [DarkSide], “Low-Mass Dark Matter Search with the DarkSide-50 Experiment,” Phys. Rev. Lett. 121, no.8, 081307 (2018) [arXiv:1802.06994 [astro-ph.HE]].
  * (52) E. Aprile et al. [XENON], “Light Dark Matter Search with Ionization Signals in XENON1T,” Phys. Rev. Lett. 123, no.25, 251801 (2019) [arXiv:1907.11485 [hep-ex]].
  * (53) M. Ackermann et al. [Fermi-LAT], “Searching for Dark Matter Annihilation from Milky Way Dwarf Spheroidal Galaxies with Six Years of Fermi Large Area Telescope Data,” Phys. Rev. Lett. 115, no.23, 231301 (2015) [arXiv:1503.02641 [astro-ph.HE]].
  * (54) A. Albert et al. [Fermi-LAT and DES], “Searching for Dark Matter Annihilation in Recently Discovered Milky Way Satellites with Fermi-LAT,” Astrophys. J. 834, no.2, 110 (2017) [arXiv:1611.03184 [astro-ph.HE]].
  * (55) M. Aaboud et al. [ATLAS], “Search for dark matter and other new phenomena in events with an energetic jet and large missing transverse momentum using the ATLAS detector,” JHEP 01, 126 (2018) [arXiv:1711.03301 [hep-ex]].
  * (56) [ATLAS], “Search for new phenomena in events with jets and missing transverse momentum in p p collisions at $\sqrt{s}$ = 13 TeV with the ATLAS detector,” ATLAS-CONF-2020-048.
  * (57) A. M. Sirunyan et al. [CMS], “Search for dark matter produced with an energetic jet or a hadronically decaying W or Z boson at $\sqrt{s}=13$ TeV,” JHEP 07, 014 (2017) [arXiv:1703.01651 [hep-ex]].
  * (58) R. Ramos, Van Que Tran and T. C. Yuan, Phys. Rev. D 103, no.7, 075021 (2021) [arXiv:2101.07115 [hep-ph]].
  * (59) P. A. Zyla et al. [Particle Data Group], “Review of Particle Physics,” PTEP 2020, no.8, 083C01 (2020).
  * (60) B. Dirgantara and C. S. Nugroho, “Effects of New Heavy Fermions on Complex Scalar Dark Matter Phenomenology in Gauged Two Higgs Doublet Model,” [arXiv:2012.13170 [hep-ph]].
  * (61) D. Feldman, Z. Liu and P. Nath, “Probing a very narrow Z-prime boson with CDF and D0 data,” Phys. Rev. Lett. 97, 021801 (2006) [hep-ph/0603039].
  * (62) M. Fabbrichesi, E. Gabrielli and G. Lanfranchi, “The Dark Photon,” [arXiv:2005.01515 [hep-ph]].
  * (63) R. Aaij et al. [LHCb], “Search for $A^{\prime}\rightarrow\mu^{+}\mu^{-}$ Decays,” Phys. Rev. Lett. 124, no.4, 041801 (2020) [arXiv:1910.06926 [hep-ex]].
  * (64) J. P. Lees et al. [BaBar], “Search for a Dark Photon in $e^{+}e^{-}$ Collisions at BaBar,” Phys. Rev. Lett. 113, no.20, 201801 (2014) [arXiv:1406.2980 [hep-ex]].
  * (65) J. R. Batley et al. [NA48/2], “Search for the dark photon in $\pi^{0}$ decays,” Phys. Lett. B 746, 178-185 (2015) [arXiv:1504.00607 [hep-ex]].
  * (66) D. Banerjee et al. [NA64], “Search for a Hypothetical 16.7 MeV Gauge Boson and Dark Photons in the NA64 Experiment at CERN,” Phys. Rev. Lett. 120, no.23, 231802 (2018) [arXiv:1803.07748 [hep-ex]].
  * (67) 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. “A Search for Short Lived Axions in an Electron Beam Dump Experiment,” Phys. Rev. Lett. 59, 755 (1987)
  * (68) J. Blümlein 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]].
  * (69) 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]].
  * (70) J. B. Dent, F. Ferrer and L. M. Krauss, “Constraints on Light Hidden Sector Gauge Bosons from Supernova Cooling,” [arXiv:1201.2683 [astro-ph.CO]];
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]];
H. K. Dreiner, J. F. Fortin, C. Hanhart and L. Ubaldi, “Supernova constraints
on MeV dark sectors from $e^{+}e^{-}$ annihilations,” Phys. Rev. D 89, no.10,
105015 (2014) [arXiv:1310.3826 [hep-ph]];
H. Merkel et al., “Search at the Mainz Microtron for Light Massive Gauge
Bosons Relevant for the Muon g-2 Anomaly,” Phys. Rev. Lett. 112, no.22, 221802
(2014) [arXiv:1404.5502 [hep-ex]];
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]];
A. Anastasi et al. [KLOE-2], “Limit on the production of a new vector boson in
$\mathrm{e^{+}e^{-}}\rightarrow{\rm U}\gamma$, U$\rightarrow\pi^{+}\pi^{-}$
with the KLOE experiment,” Phys. Lett. B 757, 356-361 (2016) [arXiv:1603.06086
[hep-ex]];
L. Marsicano, M. Battaglieri, M. Bondí, 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]];
CMS Collaboration, “Search for a narrow resonance decaying to a pair of muons
in proton-proton collisions at 13 TeV,” CMS-PAS-EXO-19-018.

  * (71) G. Aad et al. [ATLAS], “Combined measurements of Higgs boson production and decay using up to 80 fb-1 of proton-proton collision data at $\sqrt{s}=13$ TeV collected with the ATLAS experiment,” Phys. Rev. D 101, no.1, 012002 (2020) [arXiv:1909.02845 [hep-ex]].
  * (72) A. M. Sirunyan et al. [CMS], “Combined measurements of Higgs boson couplings in proton-proton collisions at $\sqrt{s}=13\,\text{TeV}$,” Eur. Phys. J. C 79, no.5, 421 (2019) [arXiv:1809.10733 [hep-ex]].
  * (73) [ATLAS], “Search for invisible Higgs boson decays with vector boson fusion signatures with the ATLAS detector using an integrated luminosity of 139 fb-1,” ATLAS-CONF-2020-008.
  * (74) G. Arcadi, M. Dutra, P. Ghosh, M. Lindner, Y. Mambrini, M. Pierre, S. Profumo and F. S. Queiroz, “The waning of the WIMP? A review of models, searches, and constraints,” Eur. Phys. J. C 78, no.3, 203 (2018) [arXiv:1703.07364 [hep-ph]].
  * (75) M. Escudero, A. Berlin, D. Hooper and M. X. Lin, “Toward (Finally!) Ruling Out Z and Higgs Mediated Dark Matter Models,” JCAP 12, 029 (2016) [arXiv:1609.09079 [hep-ph]].
  * (76) 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]].
  * (77) 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]].
  * (78) J. Alwall, M. Herquet, F. Maltoni, O. Mattelaer and T. Stelzer, “MadGraph 5 : Going Beyond,” JHEP 06, 128 (2011) [arXiv:1106.0522 [hep-ph]].
  * (79) B. Dumont, B. Fuks, S. Kraml, S. Bein, G. Chalons, E. Conte, S. Kulkarni, D. Sengupta and C. Wymant, “Toward a public analysis database for LHC new physics searches using MadAnalysis 5,” Eur. Phys. J. C 75, no.2, 56 (2015) [arXiv:1407.3278 [hep-ph]].
  * (80) G. Apollinari, O. Brüning, T. Nakamoto and L. Rossi, “High Luminosity Large Hadron Collider HL-LHC,” CERN Yellow Rep., no.5, 1-19 (2015) [arXiv:1705.08830 [physics.acc-ph]].
  * (81) M. Benedikt and F. Zimmermann, “Proton Colliders at the Energy Frontier,” Nucl. Instrum. Meth. A 907, 200-208 (2018) [arXiv:1803.09723 [physics.acc-ph]].
  * (82) N. Arkani-Hamed, T. Han, M. Mangano and L. T. Wang, “Physics opportunities of a 100 TeV proton–proton collider,” Phys. Rept. 652, 1-49 (2016) [arXiv:1511.06495 [hep-ph]].
  * (83) 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]].
  * (84) A. Alloul, N. D. Christensen, C. Degrande, C. Duhr and B. Fuks, “FeynRules 2.0 - A complete toolbox for tree-level phenomenology,” Comput. Phys. Commun. 185, 2250-2300 (2014) [arXiv:1310.1921 [hep-ph]].
  * (85) A. Belyaev, N. D. Christensen and A. Pukhov, “CalcHEP 3.4 for collider physics within and beyond the Standard Model,” Comput. Phys. Commun. 184, 1729-1769 (2013) [arXiv:1207.6082 [hep-ph]].
  * (86) 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]].
  * (87) M. Battaglieri, A. Belloni, A. Chou, P. Cushman, B. Echenard, R. Essig, J. Estrada, J. L. Feng, B. Flaugher and P. J. Fox, et al. “US Cosmic Visions: New Ideas in Dark Matter 2017: Community Report,” [arXiv:1707.04591 [hep-ph]].
  * (88) R. Agnese et al. [SuperCDMS], “Projected Sensitivity of the SuperCDMS SNOLAB experiment,” Phys. Rev. D 95, no.8, 082002 (2017) [arXiv:1610.00006 [physics.ins-det]].
  * (89) H. Ma et al. [CDEX], “CDEX Dark Matter Experiment: Status and Prospects,” J. Phys. Conf. Ser. 1342, no.1, 012067 (2020) [arXiv:1712.06046 [hep-ex]].
  * (90) E. Kou et al. [Belle-II], “The Belle II Physics Book,” PTEP 2019, no.12, 123C01 (2019) [erratum: PTEP 2020, no.2, 029201 (2020)] [arXiv:1808.10567 [hep-ex]].
  * (91) D. Banerjee et al., “Addendum to the NA64 Proposal: search for $A^{\prime}\rightarrow$ invisible and $X\rightarrow e^{+}e^{-}$ decays in 2021,” CERN-SPSC-2018-004 (SPSC-P-348-ADD-2).
  * (92) A. Caldwell, J. Chappell, P. Crivelli, E. Depero, J. Gall, S. Gninenko, E. Gschwendtner, A. Hartin, F. Keeble and J. Osborne, et al. “Particle physics applications of the AWAKE acceleration scheme,” [arXiv:1812.11164 [physics.acc-ph]].
  * (93) M. D’Onofrio, O. Fischer and Z. S. Wang, “Searching for Dark Photons at the LHeC and FCC-he,” Phys. Rev. D 101, no.1, 015020 (2020) [arXiv:1909.02312 [hep-ph]].