#  Consistency of Gauged Two Higgs Doublet Model:
Gauge Sector

Cheng-Tse Huang1, Raymundo Ramos2,
Van Que Tran2,3, Yue-Lin Sming Tsai2,4 and Tzu-Chiang Yuan2 1Interdisciplinary
Program of Sciences, National Tsing Hua University, Hsinchu 30013, Taiwan
2Institute of Physics, Academia Sinica, Nangang, Taipei 11529, Taiwan
3School of Physics, Nanjing University, Nanjing 210093, China
4Key Laboratory of Dark Matter and Space Astronomy, Purple Mountain
Observatory, Chinese Academy of Sciences, Nanjing 210008, China

###### Abstract

We study the constraints on the new parameters in the gauge sector of gauged
two Higgs doublet model using the electroweak precision test data collected
from the Large Electron Positron Collider (LEP) at and off the $Z$-pole as
well as the current Drell-Yan and high-mass dilepton resonance data from the
Large Hadron Collider (LHC). Impacts on the new parameters by the projected
sensitivities of various electroweak observables at the Circular Electron
Positron Collider (CEPC) proposed to be built in China are also discussed. We
also clarify why the Stüeckelberg mass $M_{Y}$ for the hypercharge $U(1)_{Y}$
is set to be zero in the model by showing that it would otherwise lead to the
violation of the standard charge assignments for the elementary quarks and
leptons when they couple to the massless photon.

##  I Introduction

The discovery of the 125 GeV scalar boson identified as the Higgs boson in the
Standard Model (SM) Salam:1959zz ; Glashow:1961tr ; Weinberg:1967tq ; Salam
suggested that the simple Higgs mechanism Englert:1964et ; Higgs:1964pj ;
Guralnik:1964eu  for electroweak symmetry breaking proposed by Weinberg
Weinberg:1967tq  and Salam Salam  is the choice by nature. Both Run I and Run
II data collected by the two experimental groups ATLAS and CMS at the Large
Hadron Collider (LHC) reveal no significant deviations from the SM
predictions. Alternative models for electroweak symmetry breaking like
technicolor or composite Higgs models are arguably more elegant but
necessarily more complicated. Simplicity seems to be more superior over other
criterion like complexity or elegance for model buildings.

Nevertheless experimental observations of neutrino oscillations imply there
must be new physics beyond the SM to account for the minuscule masses of
neutrinos. Missing mass problem and cosmic acceleration of our universe also
suggested the introduction of dark matter (DM) Lin:2019uvt  and dark energy
Peebles:2002gy . The standard $\Lambda$CDM model of cosmology Dodelson:2003ft
consists of the SM of particle physics plus two new ingredients, namely the
cold dark matter, which can be the weakly interacting massive particle
predicted by many new particle physics models, and a tiny positive
cosmological constant at the present time in the Einstein’s field equation for
gravity, which can be mimicked by numerous models of dark energy. Many models
of dark matter and neutrino masses require extension not only of the simple
Higgs sector but sometimes also the electroweak gauge sector of the SM as
well. Moreover, models of dark energy are often represented by new scalar
field with equation of state that can provide negative pressure in order to
explain the cosmic acceleration at late times.

Thus extension of the SM in one way or the other seems necessary if one wants
to solve the above puzzles in the neutrino sector and in cosmology. At the
same time, one should be open-minded that there might be other approaches
other than particle physics to answer some of these questions and remembering
that nature is the ultimate arbiter of all theoretical imaginations.

The gauged two Higgs doublet model (G2HDM) proposed in Huang:2015wts  was
motivated partly by the inert Higgs doublet model (IHDM) Deshpande:1977rw ;
Ma:2006km ; Barbieri:2006dq ; LopezHonorez:2006gr  of dark matter. IHDM is a
variant of the general 2HDM Branco:2011iw  with an imposed discrete $Z_{2}$
symmetry on the scalar potential and the Yukawa couplings such that one of the
Higgs doublets is odd and become a scalar dark matter candidate. Dangerous
tree level flavor changing neutral current (FCNC) interactions in the Yukawa
couplings, generally presence in the general 2HDM, are also eliminated by this
discrete symmetry. Due to its relatively simple extension of the SM, many
detailed analysis of IHDM had been done in the literature Arhrib:2013ela ;
Arhrib:2014pva ; Ilnicka:2015jba ; Belyaev:2016lok ; Eiteneuer:2017hoh ;
Borah:2018rca ; Kephart:2015oaa ; Goudelis:2013uca ; Swiezewska:2012eh ;
Arhrib:2012ia . In G2HDM, the discrete $Z_{2}$ symmetry in IHDM was not
enforced. Instead the two Higgs doublets $H_{1}$ and $H_{2}$ are grouped into
a two-dimensional irreducible representation $H=(H_{1},\,H_{2})^{\rm T}$ of a
new gauge group $SU(2)_{H}$. A priori there is no need to impose the discrete
$Z_{2}$ symmetry in G2HDM. Once we write down all renormalizable interactions
for G2HDM, this discrete symmetry emerges as an accidental symmetry
automatically. Tree level flavor changing neutral current (FCNC) interaction
in the Higgs-Yukawa couplings are also absence naturally for the SM fermions.
As long as one does not break this symmetry spontaneously, which might lead to
the domain wall problem in early universe, the $H_{2}$ doublet is naturally an
inert Higgs doublet and can play some role in dark matter physics. It is more
satisfactory to have a global discrete symmetry like the $Z_{2}$ parity that
guarantees the stability of dark matter embedded into a local symmetry. Indeed
there exists theoretical arguments showing that global continuous or discrete
symmetries are not compatible with quantum gravity Krauss:1988zc ;
Kallosh:1995hi . Detailed analysis of the complex scalar dark matter physics
in G2HDM will be presented in a forthcoming paper DMPaper .

The construction of G2HDM in Huang:2015wts  involves extension of both the
Higgs and gauge sector of the SM which we will discuss shortly in the next
section. Several phenomenological implications of G2HDM had been explored in
Huang:2015rkj ; Huang:2017bto ; Arhrib:2018sbz ; Chen:2018wjl . In particular,
we have studied recently in details the theoretical and phenomenological
constraints on the scalar sector Arhrib:2018sbz .

We note that the 2HDM augmented with an extra local abelian $U(1)_{X}$ has
been discussed in the literature Ko:2012hd ; Campos:2017dgc ; Camargo:2018klg
; Camargo:2018uzw ; Camargo:2019ukv ; Cogollo:2019mbd  to address neutrino
masses, dark matter and to avoid FCNC interactions at the tree level.

As mentioned before, all experimental data are in line with SM predictions.
The extended gauge sector of G2HDM must be challenged by electroweak precision
test (EWPT) data obtained previously at LEP-I and LEP-II as well as current
data at the LHC. Constraints must be imposed on the new parameters in the
extended gauge sector of G2HDM. The main purpose of this work is to study
these constraints on the gauge sector systematically in analogous to previous
analysis Arhrib:2018sbz  done for the scalar sector. It is also interesting to
address the sensitivities of these new parameters at the future colliders.

The contents of this paper is organized as follows: In the next Sec. II, we
review the G2HDM and highlight some of its crucial features of the gauge
sector relevant most to this work. Sec. III discusses the experimental
constraints, including the electroweak precision test constraints at and off
the $Z$-pole at LEP, Drell-Yan data from on-shell decay of the $Z$ boson at
the LHC, and the full LHC Run II data from the high-mass dilepton resonance of
an extra neutral gauge boson $Z^{\prime}$. Sec. IV contains our numerical
results from the profile likelihood analysis. We also study future
sensitivities of the new parameters in future experiments, in particular for
the Circular Electron Positron Collider (CEPC) CEPC-SPPCStudyGroup:2015csa
proposed/debated to be built in China. Finally, we summarize and conclude in
Sec. V. In Appendix A we present the formulas for the mixing angles among the
three massive neutral gauge bosons in G2HDM in terms of the fundamental
parameters in the Lagrangian of the model. The dominant two-body decay widths
for the two new neutral gauge bosons are discussed in Appendix B.

##  II G2HDM Set Up

In this section, we will start with a brief review for the set-up of G2HDM
Huang:2015wts  by specifying its particle content (Sec. II.1) and then write
down the mass spectrum of the neutral gauge bosons (Sec. II.2) and their
interactions with the SM fermions (Sec. II.3) in the model. Along the way, we
will discuss some peculiar effects for nonzero Stüeckelberg mass $M_{Y}$
associated with the hypercharge $U(1)_{Y}$.

Fields | Spin | $SU(3)_{C}$ | $SU(2)_{L}$ | $SU(2)_{H}$ | $U(1)_{Y}$ | $U(1)_{X}$
---|---|---|---|---|---|---
$H=\left(H_{1}\;\;H_{2}\right)^{\rm T}$ | 0 | 1 | 2 | 2 | $\frac{1}{2}$ | $1$
$\Delta_{H}=\left(\begin{array}[]{cc}\Delta_{3}/2&\Delta_{p}/\sqrt{2}\\\ \Delta_{m}/\sqrt{2}&-\Delta_{3}/2\end{array}\right)$ | 0 | 1 | 1 | 3 | 0 | 0
$\Phi_{H}=\left(\Phi_{1}\;\;\Phi_{2}\right)^{\rm T}$ | 0 | 1 | 1 | 2 | 0 | $1$
$Q_{L}=\left(u_{L}\;\;d_{L}\right)^{\rm T}$ | $\frac{1}{2}$ | 3 | 2 | 1 | $\frac{1}{6}$ | 0
$U_{R}=\left(u_{R}\;\;u^{H}_{R}\right)^{\rm T}$ | $\frac{1}{2}$ | 3 | 1 | 2 | $\frac{2}{3}$ | $1$
$D_{R}=\left(d^{H}_{R}\;\;d_{R}\right)^{\rm T}$ | $\frac{1}{2}$ | 3 | 1 | 2 | $-\frac{1}{3}$ | $-1$
$u_{L}^{H}$ | $\frac{1}{2}$ | 3 | 1 | 1 | $\frac{2}{3}$ | 0
$d_{L}^{H}$ | $\frac{1}{2}$ | 3 | 1 | 1 | $-\frac{1}{3}$ | 0
$L_{L}=\left(\nu_{L}\;\;e_{L}\right)^{\rm T}$ | $\frac{1}{2}$ | 1 | 2 | 1 | $-\frac{1}{2}$ | 0
$N_{R}=\left(\nu_{R}\;\;\nu^{H}_{R}\right)^{\rm T}$ | $\frac{1}{2}$ | 1 | 1 | 2 | 0 | $1$
$E_{R}=\left(e^{H}_{R}\;\;e_{R}\right)^{\rm T}$ | $\frac{1}{2}$ | 1 | 1 | 2 | $-1$ | $-1$
$\nu_{L}^{H}$ | $\frac{1}{2}$ | 1 | 1 | 1 | 0 | 0
$e_{L}^{H}$ | $\frac{1}{2}$ | 1 | 1 | 1 | $-1$ | 0
$g^{a}_{\mu}(a=1,\cdots,8)$ | 1 | 8 | 1 | 1 | 0 | 0
$W^{i}_{\mu}(i=1,2,3)$ | 1 | 1 | 3 | 1 | 0 | 0
$W^{\prime i}_{\mu}(i=1,2,3)$ | 1 | 1 | 1 | 3 | 0 | 0
$B_{\mu}$ | 1 | 1 | 1 | 1 | 0 | 0
$X_{\mu}$ | 1 | 1 | 1 | 1 | 0 | 0
Table 1: Particle content and their quantum number assignments in G2HDM.

###  II.1 Particle Content

The particle content of G2HDM is listed in Table 1 111
$u^{H}_{L},d^{H}_{L},\nu^{H}_{L},e^{H}_{L}$ in the table were denoted as
$\chi_{u},\chi_{d},\chi_{\nu},\chi_{e}$ respectively in Huang:2015wts ..
Besides the two Higgs doublets $H_{1}$ and $H_{2}$ combining to form
$H=(H_{1},H_{2})^{\rm T}$ in the fundamental representation of an extra
$SU(2)_{H}$, we introduced a triplet $\Delta_{H}$ and a doublet $\Phi_{H}$ of
this new gauge group. However $\Delta_{H}$ and $\Phi_{H}$ are singlets under
the electroweak SM gauge group $SU(2)_{L}\times U(1)_{Y}$. Only $H$ carries
both quantum numbers of the $SU(2)_{L}$ and $SU(2)_{H}$.

There are different ways of introducing new heavy fermions in the model but we
choose a simple realization: the heavy fermions together with the SM right-
handed fermions comprise $SU(2)_{H}$ doublets, while the SM left-handed
doublets are singlets under $SU(2)_{H}$. We note that heavy right-handed
neutrinos paired up with a mirror charged leptons forming $SU(2)_{L}$ doublets
was suggested before in the mirror fermion model Hung:2006ap . To render the
model anomaly-free, four additional chiral (left-handed) fermions for each
generation, all singlets under both $SU(2)_{L}$ and $SU(2)_{H}$, are included.
For the Yukawa interactions that couple among the fermions and scalars in
G2HDM, we refer our readers to Huang:2015wts  for more details, since they are
not relevant to this work.

To avoid some unwanted pieces in the scalar potential and Yukawa couplings, we
require the matter fields to carry extra local $U(1)_{X}$ charges. Thus the
complete gauge groups in G2HDM consist of $SU(3)_{C}\times SU(2)_{L}\times
U(1)_{Y}\times SU(2)_{H}\times U(1)_{X}$. Apart from the matter content of
G2HDM, there also exist the gauge bosons corresponding to the SM and the extra
gauge groups.

The salient features of G2HDM are: (i) it is free of gauge and gravitational
anomalies; (ii) renormalizable; (iii) without resorting to an ad-hoc $Z_{2}$
symmetry, an inert Higgs doublet $H_{2}$ can be naturally realized, providing
a DM candidate; (iv) due to the nonabelian $SU(2)_{H}\times U(1)_{X}$ gauge
symmetry, dangerous FCNC interactions are absent at tree level for the SM
sector; (v) the VEV of the triplet can trigger $SU(2)_{L}$ symmetry breaking
while that of $\Phi_{H}$ provides a mass to the new fermions through
$SU(2)_{H}$-invariant Yukawa couplings; etc.

###  II.2 Neutral Gauge Boson Masses

Consider the interaction basis $\\{B,W^{3},W^{\prime 3},X\\}$ for the neutral
gauge bosons and denote their mass eigenstates as $\\{A,Z_{1},Z_{2},Z_{3}\\}$.
After spontaneous symmetry breaking, the 4$\times$4 mass matrix in the
interaction basis of $\\{B,W^{3},W^{\prime 3},X\\}$ is given by Huang:2015wts

$\mathcal{M}_{\text{gauge}}^{2}=\begin{pmatrix}\frac{g^{\prime
2}v^{2}}{4}+M_{Y}^{2}&-\frac{g^{\prime}g\,v^{2}}{4}&\frac{g^{\prime}g_{H}v^{2}}{4}&\frac{g^{\prime}g_{X}v^{2}}{2}+M_{X}M_{Y}\\\
-\frac{g^{\prime}g\,v^{2}}{4}&\frac{g^{2}v^{2}}{4}&-\frac{gg_{H}v^{2}}{4}&-\frac{gg_{X}v^{2}}{2}\\\
\frac{g^{\prime}g_{H}v^{2}}{4}&-\frac{gg_{H}v^{2}}{4}&\frac{g^{2}_{H}\left(v^{2}+v^{2}_{\Phi}\right)}{4}&\frac{g_{H}g_{X}\left(v^{2}-v^{2}_{\Phi}\right)}{2}\\\
\frac{g^{\prime}g_{X}v^{2}}{2}+M_{X}M_{Y}&-\frac{gg_{X}v^{2}}{2}&\frac{g_{H}g_{X}\left(v^{2}-v^{2}_{\Phi}\right)}{2}&g_{X}^{2}\left(v^{2}+v^{2}_{\Phi}\right)+M_{X}^{2}\end{pmatrix}\;.$
(1)

Here $g$, $g^{\prime}$, $g_{H}$ and $g_{X}$ denote the gauge couplings of
$SU(2)_{L}$, $U(1)_{Y}$, $SU(2)_{H}$ and $U(1)_{X}$ respectively; $v$ and
$v_{\Phi}$ are the vacuum expectation values (VEVs) of $H_{1}$ and $\Phi_{H}$
respectively; $M_{X}$ and $M_{Y}$ are the Stüeckelberg masses for the two
abelian $U(1)_{X}$ and $U(1)_{Y}$ respectively. We note that $v_{\Delta}$ the
VEV of the triplet $\Delta_{H}$ does not enter into the neutral gauge boson
mass matrix. This is unlike the case of scalar boson mass matrix analyzed in
Arhrib:2018sbz  which involves all three VEVs, $v$, $v_{H}$ and $v_{\Delta}$.
The matrix $\mathcal{M}^{2}_{\rm gauge}$ in Eq. (1) is real and symmetric and
thus can be diagonalized by a 4$\times$4 orthogonal rotation matrix that we
will denote as ${\cal O}^{4\times 4}$

$({\cal O}^{4\times 4})^{\rm T}\cdot\mathcal{M}^{2}_{\text{gauge}}\cdot{\cal
O}^{4\times 4}=\text{diag}(0,M^{2}_{Z_{1}},M^{2}_{Z_{2}},M^{2}_{Z_{3}})\,,$
(2)

where $M^{2}_{Z_{1}}<M^{2}_{Z_{2}}<M^{2}_{Z_{3}}$. The zero mass state is
naturally identified as the photon.

Some comments on the Stüeckelberg masses $M_{X}$ and $M_{Y}$ are in order
here. It has been demonstrated in Ruegg:2003ps  that for the extension of SM
with a Stüeckelberg mass $M_{Y}$ for the hypercharge $U(1)_{Y}$, there exists
a plethora of new physical effects. Notably, besides the photon obtaining a
mass, neutrinos will couple to the photon and charged leptons will have axial
vector couplings with the photon. Nevertheless, the Stüeckelberg extension of
the SM doesn’t spoilt renormalizability of the model. All these new effects
are proportional to $M_{Y}$. Experimentally, the photon mass upper bound
deduced from modeling the solar wind in magnetohydrodynamics is
$m_{\gamma}<1\times 10^{-18}$ eV Tanabashi:2018oca , which implies $M_{Y}$
must be very tiny too. If individual Stüeckelberg mechanism is introduced for
each of the two $U(1)$s factors in G2HDM, the photon will in general obtain
nonzero mass and many results obtained in Ruegg:2003ps  apply as well. In
Huang:2015wts , we followed Kors:2005uz ; Kors:2004iz ; Kors:2004ri ;
Kors:2004dx  in which only one Stüeckelberg field was introduced for the two
factors of $U(1)$s to implement the Stüeckelberg mechanism. The matrix
$\mathcal{M}_{\text{gauge}}^{2}$ thus obtained given in Eq. (1) has zero
determinant and a massless photon can always be realized for arbitrary values
of the Stüeckelberg masses $M_{X}$ and $M_{Y}$.

In the next subsection, we will show that with a nonzero $M_{Y}$ the electric
charge assignments of the SM fermions and their heavy partners in G2HDM will
no longer be standard but instead receive milli-charge corrections like those
discussed in Ruegg:2003ps . In particular, neutrinos will couple to the photon
and all fermions also have axial vector couplings with the photon at tree
level. These peculiar effects depend on $M_{Y}$ through the mixing matrix
elements and hence necessarily small. Thus, we have strong theoretical
motivation to set $M_{Y}=0$ in what follows to avoid these unpleasant
features. For an analysis with both $M_{X}$ and $M_{Y}$ nonzero in a
Stüeckelberg $U(1)_{X}$ extension of the SM that maintains the standard QED
interaction for the SM fermions, see Feldman:2007nf ; Feldman:2007wj ;
Feldman:2006wb . The main reason why the photon-fermion couplings in G2HDM are
in general different from these previous works is due to the presence of the
extra gauge group $SU(2)_{H}$ whereas there is only one extra abelian group
$U(1)_{X}$ in Feldman:2007nf ; Feldman:2007wj ; Feldman:2006wb .

Setting $M_{Y}=0$ in G2HDM will simplify $\mathcal{M}_{\text{gauge}}^{2}$ and
allows us to write the rotation matrix in the following product form

${\cal O}^{4\times 4}_{M_{Y}=0}=\begin{pmatrix}c_{W}&-s_{W}&0&0\\\
s_{W}&c_{W}&0&0\\\ 0&0&1&0\\\
0&0&0&1\end{pmatrix}\cdot\begin{pmatrix}1&0&0&0\\\ 0&&&\\\ 0&&{\cal O}&\\\
0&&&\end{pmatrix},$ (3)

where $c_{W}$ and $s_{W}$ represent $\cos\theta_{W}$ and $\sin\theta_{W}$
respectively, with $\theta_{W}$ being the Weinberg angle defined by

$e^{\rm SM}\equiv
g\sin\theta_{W}=g^{\prime}\cos\theta_{W}=\frac{gg^{\prime}}{\sqrt{g^{2}+g^{\prime
2}}}\;.$ (4)

It is obvious that the matrix ${\cal O}^{4\times 4}_{M_{Y}=0}$ in Eq. (3) is
just the product of the SM gauge rotation matrix made into a $4\times 4$
matrix, called ${\cal O}^{4\times 4}_{\text{SM}}$, times a general $3\times 3$
orthogonal rotation matrix ${\cal O}$ which was also converted to a $4\times
4$ matrix. After applying the rotation ${\cal O}^{4\times 4}_{\text{SM}}$ to
$\mathcal{M}^{2}_{\text{gauge}}(M_{Y}=0)$, the result is

$\displaystyle{\cal O}^{4\times 4\,\rm
T}_{\text{SM}}\cdot\mathcal{M}^{2}_{\text{gauge}}(M_{Y}=0)\cdot{\cal
O}^{4\times 4}_{\text{SM}}$ $\displaystyle=\begin{pmatrix}0&0&0&0\\\
0&M_{Z^{\text{SM}}}^{2}&-\frac{g_{H}v}{2}M_{Z^{\text{SM}}}&-g_{X}vM_{Z^{\text{SM}}}\\\
0&-\frac{g_{H}v}{2}M_{Z^{\text{SM}}}&\frac{g_{H}^{2}\left(v^{2}+v_{\Phi}^{2}\right)}{4}&\frac{g_{X}g_{H}\left(v^{2}-v_{\Phi}^{2}\right)}{2}\\\
0&-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}\,,$
(5)

where $M_{Z^{\text{SM}}}=v\sqrt{g^{2}+g^{\prime 2}}/2$ is the mass of the $Z$
boson in the SM. We can consider the vanishing (1,1) element to be the mass of
the photon eigenstate $A_{\mu}$. Furthermore, according to Eqs. (2) and (3),
the remaining 3$\times$3 matrix formed by the non-vanishing elements above is
diagonalized by the orthogonal matrix ${\cal O}$. In particular, one can
parametrize ${\cal O}$ in terms of the following Tait-Bryan representation

${\cal
O}=\begin{pmatrix}c_{\psi}c_{\phi}-s_{\theta}s_{\phi}s_{\psi}&-c_{\theta}s_{\phi}&s_{\psi}c_{\phi}+s_{\theta}s_{\phi}c_{\psi}\\\
c_{\psi}s_{\phi}+s_{\theta}c_{\phi}s_{\psi}&c_{\theta}c_{\phi}&s_{\psi}s_{\phi}-s_{\theta}c_{\phi}c_{\psi}\\\
-c_{\theta}s_{\psi}&s_{\theta}&c_{\theta}c_{\psi}\end{pmatrix},$ (6)

where $s_{x}$ and $c_{x}$ stand for sine and cosine with the rotation angle
$x=\phi,\theta,\psi$ respectively. As shown in Appendix A, these rotation
angles can be represented as

$\tan({\phi})=\frac{-g_{H}vM_{Z^{\text{SM}}}(M_{X}^{2}-M_{Z_{2}}^{2}+2g_{X}^{2}v_{\Phi}^{2})}{2\Big{(}M_{Z_{2}}^{4}-\big{(}M_{Z^{\text{SM}}}^{2}+M_{X}^{2}+(v^{2}+v_{\Phi}^{2})g_{X}^{2}\big{)}M_{Z_{2}}^{2}+M_{Z^{\text{SM}}}^{2}(M_{X}^{2}+g_{X}^{2}v_{\Phi}^{2})\Big{)}}\,,$
(7)
$\tan({\theta})=\frac{-g_{X}(M_{Z_{2}}^{2}(v^{2}-v_{\Phi}^{2})+M_{Z^{\text{SM}}}^{2}v_{\Phi}^{2})}{vM_{Z^{\text{SM}}}(M_{X}^{2}-M_{Z_{2}}^{2}+2g_{X}^{2}v_{\Phi}^{2})}\sin{\phi}\,,$
(8)
$\cot({\psi})=\frac{g_{H}(M_{Z_{1}}^{2}-M_{X}^{2}-2g_{X}^{2}v_{\Phi}^{2})}{g_{X}(g_{H}^{2}v_{\Phi}^{2}-2M_{Z_{1}}^{2})}\frac{\cos{\theta}}{\sin{\phi}}-\sin{\theta}\cot{\phi}\,.$
(9)

It is easy to see that taking the limits of $g_{H}$ and $g_{X}$ go to 0, the
non-vanishing 3$\times$3 block matrix in Eq. (5) becomes ${\rm Diag}(M_{Z^{\rm
SM}}^{2},0,M_{X}^{2})$. Thus the rotation matrix $\mathcal{O}$ must be
identity. This can be realized by setting $\phi$, $\theta$ and $\psi$ to be
zeros which can be derived from Eqs. (7), (8) and (9).

We note that if one sets $M_{X}$ to zero, the mass matrix in the right-handed
side of Eq. (5) is symmetric under the interchange of $g_{H}/2\leftrightarrow
g_{X}$.

After the rotation matrix ${\cal O}$ is found, the $Z_{i}$ mass eigenstates
where $i$ runs from 1 to 3 are given by

$(Z_{1},Z_{2},Z_{3})^{\rm T}={\cal O}^{\rm T}\cdot(Z^{\rm SM},W^{\prime
3},X)^{\rm T}\,.$ (10)

The composition $Z^{\rm SM}$, $W^{\prime 3}$ and $X$ of the $Z_{i}$ mass
eigenstate is given by ${\cal O}_{1i}^{2}$, ${\cal O}_{2i}^{2}$, and ${\cal
O}_{3i}^{2}$, respectively. In general, the $Z$-pole can be any one of the
$Z_{i}$ depending on which one is actually closer to the pole by the
underlying parameter choices in G2HDM. In our analysis, we will consider there
is always at least one extra neutral gauge boson heavier than the $Z$-pole.

###  II.3 Neutral Gauge Current Interactions

The part of the Lagrangian that contains the interaction of the $Z_{i}$ with
visible matter in G2HDM is

$\mathcal{L}_{N}=g_{M}\sum_{f}\sum^{3}_{i=1}\bar{f}\gamma_{\mu}\left[\left(v_{f}^{(i)}-\gamma_{5}a^{(i)}_{f}\right)Z_{i}^{\mu}\right]f\,,$
(11)

where $g_{M}=\sqrt{g^{2}+g^{\prime 2}}/2$. The $v^{(i)}_{f}$ and $a^{(i)}_{f}$
factors are given by ($M_{Y}\neq 0$)

$\displaystyle v^{(i)}_{f}={}$ $\displaystyle\left(c_{W}{\cal O}^{4\times
4}_{2,i+1}-s_{W}{\cal O}^{4\times 4}_{1,i+1}\right)T^{3}_{f}+2Q_{f}s_{W}{\cal
O}^{4\times 4}_{1,i+1}$ $\displaystyle+\frac{1}{\sqrt{g^{2}+g^{\prime
2}}}\left(X_{R}g_{X}{\cal O}^{4\times 4}_{4,i+1}+T^{3H}_{f_{R}}g_{H}{\cal
O}^{4\times 4}_{3,i+1}\right),$ (12) $\displaystyle a^{(i)}_{f}={}$
$\displaystyle\left(c_{W}{\cal O}^{4\times 4}_{2,i+1}-s_{W}{\cal O}^{4\times
4}_{1,i+1}\right)T^{3}_{f}$ $\displaystyle-\frac{1}{\sqrt{g^{2}+g^{\prime
2}}}\left(X_{R}g_{X}{\cal O}^{4\times 4}_{4,i+1}+T^{3H}_{f_{R}}g_{H}{\cal
O}^{4\times 4}_{3,i+1}\right).$ (13)

Here $T^{3}_{f}$ is the $SU(2)_{L}$ isospin charge and $Q_{f}$ is the electric
charge in units of $e^{\rm SM}$ for the SM fermion $f$ where $e^{\rm SM}$ is
given by Eq. (4). They are related to the $U(1)_{Y}$ hypercharge by the
standard formula $Q^{\rm SM}_{f}=T^{3}_{f}+Y_{f}$. The charges due to the new
gauge symmetries are $X_{R}$ as the $U(1)_{X}$ charge of the corresponding
$f_{R}$ and $T^{3H}_{f_{R}}$ is the $SU(2)_{H}$ analogous of the $SU(2)_{L}$
isospin $T^{3}$ again for the corresponding $f_{R}$. We simply define
$T^{3H}_{f_{R}}=\pm 1/2$ depending on $f_{R}$ belongs to the upper or lower
component of an $SU(2)_{H}$ doublet.

For the photon-fermion couplings in G2HDM, we obtain

${\cal L}_{\gamma}=-e^{\text{SM}}\sum_{f}\bar{f}\gamma_{\mu}\left(Q^{\rm
G2HDM}_{f}-a^{\gamma}_{f}\gamma_{5}\right)A^{\mu}f\;,$ (14)

where

$\displaystyle Q^{\rm G2HDM}_{f}={}$ $\displaystyle\frac{{\cal O}^{4\times
4}_{1,1}}{c_{W}}Q^{\rm SM}_{f}+\frac{T^{3}_{f}}{2}\left(\frac{{\cal
O}^{4\times 4}_{2,1}}{s_{W}}-\frac{{\cal O}^{4\times 4}_{1,1}}{c_{W}}\right)$
$\displaystyle+\frac{1}{2e^{\text{SM}}}\left(g_{X}{\cal O}^{4\times
4}_{4,1}X_{R}+g_{H}{\cal O}^{4\times 4}_{3,1}T^{3H}_{f_{R}}\right)\;,$ (15)
$\displaystyle a^{\gamma}_{f}={}$
$\displaystyle\frac{T^{3}_{f}}{2}\left(\frac{{\cal O}^{4\times
4}_{2,1}}{s_{W}}-\frac{{\cal O}^{4\times 4}_{1,1}}{c_{W}}\right)$
$\displaystyle-\frac{1}{2e^{\text{SM}}}\left(g_{X}X_{R}{\cal O}^{4\times
4}_{4,1}+g_{H}T^{3H}_{f_{R}}{\cal O}^{4\times 4}_{3,1}\right)\;.$ (16)

Thus, with both nonzero $M_{X}$ and $M_{Y}$, the electromagnetism interaction
in G2HDM is in general different from the SM case. The standard charge
assignment for every SM fermion will suffer from an overall correction factor
of $\mathcal{O}^{4\times 4}_{1,1}/c_{W}$ plus two correction terms, and there
is also a non-vanishing axial vector coupling.

Next, we can take the limit $M_{Y}=0$ and write the corresponding expressions.
By replacing the elements of ${\cal O}^{4\times 4}$ by ${\cal O}^{4\times
4}_{1,1}={\cal O}^{4\times 4}_{2,2}=c_{W}$ and $-{\cal O}^{4\times
4}_{1,2}={\cal O}^{4\times 4}_{2,1}=s_{W}$ as in Eq. (3), one can find the
following new expressions for the vector and axial vector couplings

$\displaystyle v^{(i)}_{f(M_{Y}=0)}={}$
$\displaystyle(T^{3}_{f}-2Q_{f}s_{W}^{2}){\cal
O}_{1i}+\frac{1}{\sqrt{g^{2}+g^{\prime 2}}}\left(X_{R}g_{X}{\cal
O}_{3i}+T^{3H}_{f_{R}}g_{H}{\cal O}_{2i}\right),$ (17) $\displaystyle
a^{(i)}_{f(M_{Y}=0)}={}$ $\displaystyle T^{3}_{f}\,{\cal
O}_{1i}-\frac{1}{\sqrt{g^{2}+g^{\prime 2}}}\left(X_{R}g_{X}{\cal
O}_{3i}+T^{3H}_{f_{R}}g_{H}{\cal O}_{2i}\right).$ (18)

Similarly, one can do the same substitutions on Eqs. (II.3) and (II.3)
together with ${\cal O}^{4\times 4}_{3,1}={\cal O}^{4\times 4}_{4,1}=0$ and
check that the photon coupling to the SM fermions goes back to the SM
expression $Q^{\rm G2HDM}_{f}=Q^{\rm SM}_{f}=T^{3}_{f}+Y_{f}$ while all the
axial vector couplings $a^{\gamma}_{f}$ vanish. This is the main physical
reason why we set $M_{Y}=0$ so as to reproduce the standard photon-fermion
couplings. For $M_{X}$, it can be arbitrary and is naturally to consider the
light and heavy scenarios where it is smaller and greater than the $Z$-boson
mass respectively.

Obviously, the formulas obtained in this subsection for the couplings of the
neutral gauge bosons with the SM fermions also hold for the heavy fermions in
G2HDM.

##  III The Constraints

###  III.1 Constraints from Precision Electroweak Data at LEP-I

The interaction of $Z$ boson with SM fermions is described by the Lagrangian
in Eq. (11). For the case of $M_{Y}=0$ limit, the tree-level couplings are
shown in Eqs. (17) and (18). For more precise calculation, we include the
radiation corrections from propagator self-energies and flavor specific vertex
corrections to the $Z$ boson and fermions couplings Erler:2004nh ;
ALEPH:2005ab , which now are given by 222We ignore loop corrections related to
the new gauge couplings $g_{H}$ and $g_{X}$. (suppressing $M_{Y}=0$ in the
subscripts)

$\displaystyle v_{f}^{i}={}$
$\displaystyle\sqrt{\rho_{f}}\,(T^{3}_{f}-2\kappa_{f}Q_{f}s_{W}^{2}){\cal
O}_{1i}+\frac{1}{\sqrt{g^{2}+g^{\prime 2}}}\left(X_{R}g_{X}{\cal
O}_{3i}+T^{3H}_{f_{R}}g_{H}{\cal O}_{2i}\right)\,,$ (19) $\displaystyle
a_{f}^{i}={}$ $\displaystyle\sqrt{\rho_{f}}\,T^{3}_{f}\,{\cal
O}_{1i}-\frac{1}{\sqrt{g^{2}+g^{\prime 2}}}\left(X_{R}g_{X}{\cal
O}_{3i}+T^{3H}_{f_{R}}g_{H}{\cal O}_{2i}\right)\,,$ (20)

where $i$ in this work is either equal to 1 or 2 depending which mass
eigenstate is closest to $Z$-pole. The parameters $\rho_{f}$ and $\kappa_{f}$
are loop corrections quantities. The decay of the $Z$ boson into fermions and
anti-fermions in the on-shell renormalization scheme is given by Baur:2001ze ;
Erler:2004nh

$\Gamma(Z\rightarrow
f\bar{f})=N_{f}^{c}\Gamma_{o}\mathcal{R}_{f}\sqrt{1-4\mu_{f}^{2}}\Bigg{[}|v_{f}|^{2}(1+2\mu_{f}^{2})+|a_{f}|^{2}(1-4\mu_{f}^{2})\Bigg{]}\,,$
(21)

where $N_{f}^{c}$ is the color factor (1 for leptons and 3 for quarks),
$\Gamma_{o}=G_{F}M^{3}_{Z}/6\sqrt{2}\pi$, $\mu_{f}=m_{f}/M_{Z}$ and

$\mathcal{R}_{f}=\left(1+\delta_{f}^{QED}\right)\left(1+\frac{N_{f}^{c}-1}{2}\delta_{f}^{QCD}\right)\,,$
(22)

with

$\displaystyle\delta_{f}^{QED}$ $\displaystyle=$
$\displaystyle\frac{3\alpha}{4\pi}Q_{f}^{2}\,,$ (23)
$\displaystyle\delta_{f}^{QCD}$ $\displaystyle=$
$\displaystyle\frac{{\alpha_{s}}}{\pi}+1.409\left({\frac{{\alpha_{s}}}{\pi}}\right)^{2}-12.77\left({\frac{{\alpha_{s}}}{\pi}}\right)^{3}-Q_{f}^{2}\frac{{\alpha\alpha_{s}}}{{4\pi^{2}}}\,.$
(24)

Here $Q_{f}$ is the electric charge of the fermion $f$ in unit of $e^{\rm
SM}$, and $\alpha$ and $\alpha_{s}$ are the fine-structure and strong coupling
constants, respectively, evaluated at the $M_{Z}$ scale. It is understood that
the couplings $v_{f}$ and $a_{f}$ in Eq. (21) should be replaced by
$v^{i}_{f}$ and $a^{i}_{f}$ in Eqs. (19) and (20) respectively with $i=1$ or 2
depending which $M_{Z_{i}}$ is closest to the $Z$-pole $M_{Z}$.

We also investigate some $Z$-pole ($\sqrt{s}\approx M_{Z}$) observables,
including the ratio of partial decay width of $Z$ boson

$R_{l}=\frac{\Gamma_{\rm had}}{\Gamma_{l^{+}l^{-}}}\,,\hskip
28.45274ptR_{q}=\frac{\Gamma_{q\bar{q}}}{\Gamma_{\rm had}}\,,$ (25)

the hadronic cross-section

$\sigma_{\rm had}=\frac{12\pi\Gamma_{e^{+}e^{-}}\Gamma_{\rm
had}}{M_{Z}^{2}\,\Gamma_{Z}^{2}}\,,$ (26)

the parity violation quantity

$A_{f}=\frac{2v_{f}a_{f}}{v^{2}_{f}+a^{2}_{f}}\,,$ (27)

and the forward-backward asymmetry quantity

$A_{\rm FB}=\frac{3}{4}A_{f}\frac{A_{e}+P_{e}}{1+P_{e}A_{e}}\,,$ (28)

where $P_{e}$ is the initial $e^{-}$ polarization. Recall that at LEP-I
$P_{e}=0$, in this case

$A^{(0,f)}_{\rm FB}=\frac{3}{4}A_{e}A_{f}\,.$ (29)

A summary of the electroweak observables at $Z$-pole from various experiments
Tanabashi:2018oca  is presented in Table 2.

Observables | LEP Data | CEPC Precision CEPC-SPPCStudyGroup:2015csa  | Standard Model
---|---|---|---
$M_{Z}$ [GeV] | 91.1876 $\pm$ 0.0021 | $5\times 10^{-4}$ | 91.1884 $\pm$ 0.0020
$\Gamma_{Z}$ [GeV] | 2.4952 $\pm$ 0.0023 | $5.06\times 10^{-4}$ | 2.4942 $\pm$ 0.0008
$\Gamma_{\rm had}$ [GeV] | 1.7444 $\pm$ 0.0020 | — | 1.7411 $\pm$ 0.0008
$\Gamma_{\rm inv}$ [MeV] | 499.0 $\pm$ 1.5 | — | 501.44 $\pm$ 0.04
$\Gamma_{l^{+}l^{-}}$ [MeV] | 83.984 $\pm$ 0.086 | — | 83.959 $\pm$ 0.008
$\sigma_{\rm had}[{\rm nb}]$ | 41.541 $\pm$ 0.037 | — | 41.481 $\pm$ 0.008
$R_{e}$ | 20.804 $\pm$ 0.050 | — | 20.737 $\pm$ 0.010
$R_{\mu}$ | 20.785 $\pm$ 0.033 | $0.05\%$ | 20.737 $\pm$ 0.010
$R_{\tau}$ | 20.764 $\pm$ 0.045 | $0.05\%$ | 20.782 $\pm$ 0.010
$R_{b}$ | 0.21629 $\pm$ 0.00066 | $0.08\%$ | 0.21582 $\pm$ 0.00002
$R_{c}$ | 0.1721 $\pm$ 0.0030 | — | 0.17221 $\pm$ 0.00003
$A^{(0,e)}_{\rm FB}$ | 0.0145 $\pm$ 0.0025 | — | 0.01618 $\pm$ 0.00006
$A^{(0,\mu)}_{\rm FB}$ | 0.0169 $\pm$ 0.0013 | — | 0.01618 $\pm$ 0.00006
$A^{(0,\tau)}_{\rm FB}$ | 0.0188 $\pm$ 0.0017 | — | 0.01618 $\pm$ 0.00006
$A^{(0,b)}_{\rm FB}$ | 0.0992 $\pm$ 0.0016 | $0.15\%$ | 0.1030 $\pm$ 0.0002
$A^{(0,c)}_{\rm FB}$ | 0.0707 $\pm$ 0.0035 | — | 0.0735 $\pm$ 0.0001
$A^{(0,s)}_{\rm FB}$ | 0.0976 $\pm$ 0.0114 | — | 0.1031 $\pm$ 0.0002
$A_{e}$ | 0.15138 $\pm$ 0.00216 | — | 0.1469 $\pm$ 0.0003
$A_{\mu}$ | 0.142 $\pm$ 0.015 | — | 0.1469 $\pm$ 0.0003
$A_{\tau}$ | 0.136 $\pm$ 0.015 | — | 0.1469 $\pm$ 0.0003
$A_{b}$ | 0.923 $\pm$ 0.020 | — | 0.9347
$A_{c}$ | 0.670 $\pm$ 0.027 | — | 0.6677 $\pm$ 0.0001
$A_{s}$ | 0.0895 $\pm$ 0.091 | — | 0.9356
Table 2: The electroweak observables at the $Z$-pole. The second, third and
last column are the LEP measurement Tanabashi:2018oca , CEPC preliminary
conceptual design report CEPC-SPPCStudyGroup:2015csa , and the SM prediction
Tanabashi:2018oca , respectively.

From the data in Table 2, we build the Chi-squared for the electroweak
observables at $Z$-pole as follows

$\displaystyle\chi^{2}_{Z-{\rm pole}}$ $\displaystyle=$
$\displaystyle\chi^{2}_{M_{Z}}+\chi^{2}_{\sigma_{\rm
had}}+\max\left[\chi^{2}_{\Gamma_{Z}},(\chi^{2}_{\Gamma_{\rm
had}}+\chi^{2}_{\Gamma_{\rm inv}}+\chi^{2}_{\Gamma_{l^{+}l^{-}}})\right]$ (30)
$\displaystyle+\sum_{f=(e,\mu,\tau,b,c)}{\chi^{2}_{R_{f}}}+\sum_{f=(e,\mu,\tau,b,c,s)}{(\chi^{2}_{A_{f}}+\chi^{2}_{A_{\rm
FB}^{(0,f)}})}\;.$

Note that we have considered the correlations between the total decay width of
$Z$ boson and its partial decay widths to hadrons, invisibles and dilepton.
For each $\chi^{2}_{i}$ on the right-handed side of Eq. (30), it is given by
the standard expression, namely

$\displaystyle\chi^{2}_{i}=\frac{\left({O}_{i}^{\rm exp}-O_{i}^{\rm
th}\right)^{2}}{\left(\Delta O_{i}^{\rm exp}\right)^{2}}\;,$ (31)

where $O_{i}^{\rm exp/th}$ represents the experimental/theoretical value of
any one of the 23 electroweak observables listed in Table 2 and $\Delta O^{\rm
exp}_{i}$ is the corresponding experimental uncertainty.

###  III.2 Contact Interactions at LEP-II

We also include constraints from data above the $Z$-pole by considering the
LEP-II measurements related to contact interactions taking the following form
of effective Lagrangian

$\mathcal{L}_{\mathrm{eff}}=\frac{\pm 4\pi}{(1+\delta_{ef})(\Lambda^{\pm
f}_{\alpha\beta})^{2}}\left(\bar{e}\gamma^{\mu}P_{\alpha}e\bar{f}\gamma_{\mu}P_{\beta}f\right)\;,$
(32)

where $P_{\alpha,\beta}$ represent the chirality projection operators with
$\alpha,\beta$ being $L$ or $R$ for left-handed or right-handed fermions,
respectively. The sign of Eq. (32) depends on whether the interference between
the contact interaction it parametrizes and the SM process is constructive
($+$) or destructive ($-$). There is a total of 6 combinations for the
$\alpha\beta$ indices of $\Lambda^{\pm f}_{\alpha\beta}$:
$\alpha\beta=\\{{LL,LR,RL,RR,VV,AA}\\}$, which are also called models. The
limits on $\Lambda^{\pm f}_{\alpha\beta}$ set by LEP-II are given in Table
3.15 of Ref. Schael:2013ita . The strongest constraint is given by
$\Lambda^{+l}_{VV}>24.6$ TeV. By using these $\Lambda^{\pm f}_{\alpha\beta}$
values, we are able to reconstruct the cross section for new physics processes
based on the Lagrangian in Eq. (32).

To improve the analysis of this section, in particular for the cases where the
mass of one of the gauge bosons is below the $Z$-pole, we calculate the
additional $Z$-like mediator contribution 333 In what follows, we will denote
the extra neutral gauge boson as $Z^{\prime}$ or $Z_{i}$ depending on whether
we refer to the experimental data or G2HDM. to the $e^{-}e^{+}\rightarrow
Z_{i}\rightarrow f\bar{f}$ scattering cross section. In the case $f=e$ we have
the contribution of both $s$ and $t$ channels while for $f\neq e$ only the $s$
channel contributes. Note that here we do not need the SM contributions such
as the photon and $Z$ exchange not considered in Eq. (32). In the massless
approximation for all the external fermions, the amplitudes for the $s$ and
$t$ channels and for the interference term between them are given by:

$\displaystyle|\mathcal{M}_{s}|^{2}$
$\displaystyle=\frac{2g_{M}^{4}\left\\{\left[(a_{f}^{i})^{4}+(v^{i}_{f})^{4}\right]\left(s^{2}+2st+2t^{2}\right)-2(a^{i}_{f}v^{i}_{f})^{2}\left(s^{2}+2st-2t^{2}\right)\right\\}}{\left(M_{Z_{i}}^{2}-s\right)^{2}}\,,$
(33) $\displaystyle|\mathcal{M}_{t}|^{2}$
$\displaystyle=\frac{2g_{M}^{4}\left[(a^{i}_{f})^{4}\left(s^{2}+t^{2}\right)-2(a^{i}_{f}v^{i}_{f})^{2}\left(s^{2}-3t^{2}\right)+(v^{i}_{f})^{4}\left(s^{2}+t^{2}\right)\right]}{\left(M_{Z_{i}}^{2}+s+t\right)^{2}}\,,$
(34) $\displaystyle|\mathcal{M}_{st}|^{2}$
$\displaystyle=\frac{4g_{M}^{4}t^{2}\left[(a^{i}_{f})^{4}+6(a^{i}_{f}v^{i}_{f})^{2}+(v^{i}_{f})^{4}\right]}{\left(M_{Z_{i}}^{2}-s\right)\left(M_{Z_{i}}^{2}+s+t\right)}\,,$
(35)

where $s$ is the center of mass energy squared, $t=s(\cos\varphi-1)/2$ and
$\varphi$ is the angle between incoming and outgoing particles. This angle
$\varphi$ should be integrated to obtain the final cross section. The
resulting cross section has to be compared against the cross section obtained
using the effective Lagrangian in Eq. (32) with the $\Lambda^{\pm
f}_{\alpha\beta}$ given by the experimental result. The couplings $v^{i}_{f}$
and $a^{i}_{f}$ have $i=$ 1 or 2 depending on whether we are analyzing light
or heavy $M_{X}$ scenario. For $i=3$, we assume $M_{Z_{3}}$ is much heavier
than $M_{Z}$ so that its contributions are negligible. To be able to construct
a $\chi^{2}$ from the LEP-II 95% C.L. limit, we calculate the corresponding
95% C.L. cross section and compare against the theoretical result. When our
theoretical result matches the 95% C.L. with null-signal assumption, the
corresponding $\chi^{2}$ value should be 2.71 444 For a Gaussian distribution,
the value of $\Delta\chi^{2}=2.71$ corresponds to the 90% C.L. of a two-tailed
test, but it also equivalent to the 95% C.L. of a one-tailed test that we are
using. . In this case, we calculate the $\chi^{2}$ value using

$\chi^{2}_{\rm LEP-
II}=2.71\times\left[\frac{\sigma_{\text{G2HDM}}\left(e^{+}e^{-}\rightarrow
Z_{i}\rightarrow f\bar{f}\right)}{\overline{\sum}\sigma_{\rm
eff}\left(\Lambda^{\pm f\text{(95\%)}}_{\alpha\beta}\right)}\right]^{2}\;,$
(36)

where $\sigma_{\text{eff}}$ is the cross section obtained using the effective
Lagrangian of Eq. (32) with the experimental results for $\Lambda^{\pm
f}_{\alpha\beta}$ given in Ref. Schael:2013ita  for different combinations of
the chirality. The effective cross sections for different combinations of
$\alpha\beta=\\{LL,RR,LR,RL\\}$ from the data are summed and averaged. We do
not consider the combinations of VV and AA since they are not independent from
the other polarizations considered above. Note that Eq. (36) goes to zero when
the theoretical cross section vanishes (SM limit) as one would expect.

In the light $M_{X}$ scenario (see Sec. IV.3) in which one of the new neutral
gauge boson is too light and invalidates the effective contact interaction
approach, it is mandatory to recast the LEP-II constraints for the contact
interactions into the cross section level to do the analysis. We checked that
for the heavy $M_{X}$ scenario, using either the effective contact interaction
or cross section approach give the same results.

###  III.3 Drell-Yan Constraints at the LHC

In this section we recap the experiments of the Drell-Yan cross section for SM
$Z$-boson and heavy $Z^{\prime}$ at the LHC.

####  III.3.1 $Z$-boson on-shell decay at the LHC

By using the measurement of the Drell-Yan cross section for the $Z$-boson
production, the properties of the $Z$ are well determined at the LHC. Among
all the final states of the $Z$-boson decay, the dilepton signature is the
most relevant to distinguish signal from background. It is commonly believed
that the Drell-Yan constraint $q\bar{q}\rightarrow Z\rightarrow l^{+}l^{-}$
from the LHC is weaker than LEP EWPT data because of the relative larger
uncertainties from the hadronic background than the QED background. However,
to be careful, we first check a direct Drell-Yan constraints from the LHC
Aaboud:2017buh . The data of electron-positron pair ($ee$) and muon-pair
($\mu\mu$) final states are given by Tables 3 and 4 respectively in Ref.
Aaboud:2017buh . In the signal region located around $Z$-boson mass (the
invariant mass $80<m_{ll}/\;\mathrm{GeV}<120$), we found that the systematic
uncertainties of Drell-Yan background is larger than the data statistic
uncertainties in both $ee$ or $\mu\mu$ final state. We have also checked that
the EWPT constraints in Table 2 are much stronger than LHC Drell-Yan
constraint.

(a)

(b)

Figure 1: The Feynman diagrams of electroweak Drell-Yan process.

On the other hand, $Z$-boson can be singly produced either by radiation from
the incoming partons (Fig. 1a) or $t$-channel exchange of a $W$ gauge boson
(Fig. 1b). To constrain the G2HDM modified $Zl^{+}l^{-}$ couplings, the later
process is more useful than the former because QCD processes usually suffer
from larger systematical uncertainties than the electroweak ones. Recently,
ATLAS Aaboud:2017emo  reported a fiducial electroweak cross section of
$\sigma^{Zjj}_{\rm{EW}}=119\pm 16\pm 20\pm 2$ fb and
$\sigma^{Zjj}_{\rm{EW}}=34.2\pm 5.8\pm 5.5\pm 0.7$ fb for dijet invariant
masses $m_{jj}$ greater than $250\;\mathrm{GeV}$ and $1\;\mathrm{TeV}$,
respectively. The SM simulated cross sections $\sigma^{Zjj}_{\rm{EW}}({\rm
SM})$ are also given in Table 5 of Ref. Aaboud:2017emo , where central values
and the uncertainties are given as $125.2\pm 3.4$ fb for
$m_{jj}>250\;\mathrm{GeV}$ and $38.5\pm 1.5$ fb for $m_{jj}>1\;\mathrm{TeV}$.

Comparing with the SM, except for the $Z{l^{+}l^{-}}$ couplings, the G2HDM did
not modify much of the cross section. Namely, the electroweak cross section of
the G2HDM version can be simply rescaled as

$\sigma^{Zjj}_{\rm{EW}}({\rm G2HDM})=\sigma^{Zjj}_{\rm{EW}}({\rm
SM})\times\mathcal{R}\,,$ (37)

where

$\mathcal{R}=\left[\frac{\mathcal{C}^{ZWW}_{\rm{G2HDM}}}{\mathcal{C}^{ZWW}_{\rm{SM}}}\right]^{2}\frac{{\rm
BR}^{\rm{G2HDM}}_{Z\rightarrow f\bar{f}}}{{\rm BR}^{\rm{SM}}_{Z\rightarrow
f\bar{f}}}={\cal O}_{11}^{2}\frac{{\rm BR}^{\rm{G2HDM}}_{Z\rightarrow
f\bar{f}}}{{\rm BR}^{\rm{SM}}_{Z\rightarrow f\bar{f}}}\,,$ (38)

and $f=e,\mu$. However, similar to direct Drell-Yan $Z$ boson search, we found
that the value of $\mathcal{R}$ is not easy to derivate from unity and the
power of constraining the parameter space in G2HDM is not as strong as LEP
EWPT constraints.

Finally, we have numerically verified that the allowed G2HDM parameter space
is hardly changed at all whether the direct and electroweak Drell-Yan $Z$
boson constraints at the LHC are included or not. Again, this is because both
constraints at the LHC are much weaker than LEP EWPT constraints. Hence, we
will not take into account the LHC Drell-Yan constraints from the on-shell $Z$
decay in our numerical works so as to save some computer resources.

####  III.3.2 LHC $Z^{\prime}$ Boson Search at High-mass Dilepton Resonances

The Drell-Yan constraints can also be very powerful for the new gauge bosons
in G2HDM once they can be singly produced Huang:2017bto . Unlike the study in
Ref. Huang:2017bto  which only $W^{{}^{\prime}3}$-like $Z_{i}$ is considered,
we extend it here to any $Z_{i}$ with all the possible composition. Recently,
ATLAS collaboration Aad:2019fac  reported a new result on dilepton resonances
with the integrated luminosity of $139$ fb-1 and the center-of-mass energy
$\sqrt{s}=13\;\mathrm{TeV}$. They indicated that the lower limit on the mass
of $Z^{\prime}$ boson for a simplified model can be raised up to
$4-5\;\mathrm{TeV}$. Considering this new measurement, we update the
constraints of the heavy neutral gauge boson masses in G2HDM and the upper
limits of $g_{H}$ and $g_{X}$.

In Fig. 3 of Ref. Aad:2019fac , one can see the upper limits of cross section
times branching ratio BR($Z^{\prime}\rightarrow l^{+}l^{-}$) are based on the
ratio of the total width $\Gamma_{Z^{\prime}}$ of $Z^{\prime}$ divided by its
$M_{Z^{\prime}}$. Depending on this ratio, the limits can be altered by a
factor of $\sim 5$. As shown in Appendix B, the $\Gamma_{Z_{i}}/M_{Z_{i}}$ in
the G2HDM shall be always less than $0.06$. Hence, for a conservative
approach, we can simply apply the ATLAS result by using their upper limit
associated with $\Gamma_{Z^{\prime}}/M_{Z^{\prime}}=0.06$.

Furthermore, the $Z_{i}$ total decay width relies on whether $Z_{i}$ decays to
the new particles in G2HDM. The heavy new fermions in G2HDM are assumed to be
very heavy so that they do not affect the EW-scale physics in any significant
way. On the other hand, the $Z_{i}$ invisible decay to a scalar DM pair can be
a more important channel because the upper limits of various parameters can be
weaker than the one without taking into account the $Z_{i}$ decays to the DM
pair. The openings of the scalar channels as well as other channels with one
vector and one scalar particles in the final states of $Z_{i}$ decay makes the
parameter spaces of the gauge and scalar sectors entangle with each other.
Thus a complete analysis becomes quite formidable. In Eq. (76), one can see
that the invisible decay width of $Z_{i}$ has two different limits, $M_{D}\ll
M_{Z_{i}}$ for maximum invisible decay and $M_{D}>M_{Z_{i}}$ for zero
invisible decay. For the sake of simplicity, we will be contented by
presenting the results based on these two benchmark invisible decay widths. In
this study, we adopt $M_{D}=M_{Z_{i}}/10$ for maximum invisible decay but we
found that the $\Gamma(Z_{i}\rightarrow DD^{*})$ can differ within an accepted
range of $\sim 6\%$ comparing with the massless $M_{D}$ case.

Using MadGraph5 Alwall:2014hca , we compute the cross section
$\sigma(pp\rightarrow Z_{i})$. Since we enforce that the cross section is
computed at the resonance, we only used a minimum cut given by the default
parameter card in MadGraph5. It is very CPU time consuming to estimate the
cross section point by point throughout all the parameter space. Nevertheless,
the cross section can be obtained by simply rescaling the vector and axial
vector couplings $v^{i}_{f}$ and $a^{i}_{f}$ using Eqs. (17) and (18). Hence,
by using the same reasoning as before we include the latest ATLAS $Z^{\prime}$
limit in our scan by using the following chi-squared function

$\chi^{2}_{\rm{ATLAS}}=2.71\times\left[\frac{\sigma_{\text{G2HDM}}(pp\rightarrow
Z_{i})\times{\rm BR}_{\text{G2HDM}}(Z_{i}\rightarrow
l^{+}l^{-})}{\sigma^{\rm{95\%}}_{\text{ATLAS}}(pp\rightarrow
Z^{\prime})\times{\rm BR}^{\rm{95\%}}_{\text{ATLAS}}(Z^{\prime}\rightarrow
l^{+}l^{-})}\right]^{2}\;,$ (39)

where the branching ratio ${\rm BR}_{\text{G2HDM}}(Z_{i}\rightarrow
l^{+}l^{-})$ can be found in Appendix B and
$\sigma^{\rm{95\%}}_{\text{ATLAS}}(pp\rightarrow Z^{\prime})\times{\rm
BR}^{\rm{95\%}}_{\text{ATLAS}}(Z^{\prime}\rightarrow l^{+}l^{-})$ is $95\%$
C.L. taken from the curve associated with
$\Gamma_{Z^{\prime}}/M_{Z^{\prime}}=0.06$ in Fig. 3 of Ref. Aad:2019fac .

##  IV Results

###  IV.1 Numerical Methodology

Our aim is to determine the $68\%$ and $95\%$ allowed parameter space of the
G2HDM which are favored by all of the experimental data presented in the
previous section. In this paper, we will use the profile-likelihood (PL)
method to perform the statistical data analysis. We recap the PL method in the
following. Briefly, the PL method is a well popular statistical method to deal
with the multi-dimensional parameter space which treats the unwanted
parameters as nuisance parameters. In other words, if a proposed model has
$n$-dimensional parameter space and we are only interested in $p$ of those
dimensions, then the PL method can remove the unwanted $n-p$ dimensions which
we are not interested in, by maximizing the likelihood over them.

There are 4 new parameters in the gauge sector of G2HDM. They are the two new
gauge couplings $g_{H}$ and $g_{X}$ and the two new scales $v_{\Phi}$ and
$M_{X}$. Our results will be presented in two-dimensional parameter regions
with 68% and 95% confidence levels (C.L.). Take the plane ($g_{H},g_{X}$) as
an example. After marginalizing over the other two parameters $v_{\Phi}$ and
$M_{X}$, an integral of the likelihood function $\mathcal{L}(g_{H},g_{X})$ can
be written as

$\frac{\int_{\mathcal{C}}\mathcal{L}(g_{H},g_{X})dg_{H}dg_{X}}{\rm{normalization}}=\varrho\;,$
(40)

where $\mathcal{C}$ is the smallest area bound with a fraction $\varrho$ of
the total probability and the normalization in the denominator is the total
probability with $\mathcal{C}\rightarrow\infty$.

The total $\chi^{2}_{\rm Total}(g_{H},g_{X},v_{\Phi},M_{X})$ we will use in
our analysis is the sum of Eqs. (30), (36), and (39), namely

$\chi^{2}_{\rm Total}=\chi^{2}_{Z-{\rm pole}}+\chi^{2}_{\rm LEP-
II}+\chi^{2}_{\rm ATLAS}\;,$ (41)

where we have suppressed the arguments of all the $\chi^{2}$ functions. We
adopt the statistical sensitivity as

$\Delta\chi^{2}=\chi^{2}_{{\rm Total}}-{\rm min}\,(\chi^{2}_{{\rm Total}})\;.$
(42)

Since our likelihood is modeled as a pure Gaussian distribution, i.e.
$\mathcal{L}\propto\exp(-\chi^{2}/2)$, one can connect the $\chi^{2}$ to the
confidence level: the 68% (95%) C.L. in a two dimensional parameter space
corresponding to
$\Delta\chi^{2}=-2\ln(\mathcal{L}/\mathcal{L}_{\rm{max}})=2.30~{}(5.99)$. Here
$\mathcal{L}_{\rm{max}}$ is the maximum value of the likelihood in the region
$\mathcal{C}$.

There are two interesting scenarios: (i) heavy $M_{X}$ and (ii) light $M_{X}$.
The heavy $M_{X}$ scenario will result in two new heavy neutral gauge bosons
$Z_{2}\equiv Z^{\prime}$ and $Z_{3}\equiv Z^{\prime\prime}$, and the measured
boson located at $Z$-pole will be the lightest one, $Z_{1}\equiv Z$. However,
the light $M_{X}$ scenario will result in a new boson $Z_{1}$ lighter than the
$Z$-pole which is usually called dark $Z$ ($Z_{D}$) or dark photon
($\gamma_{D}$). In this case, $Z_{2}$ corresponds to the $Z$-pole $Z_{2}\equiv
Z$ and $Z_{3}\equiv Z^{\prime}$. Hence, we choose our $M_{X}$ scan ranges for
two scenarios,

$\frac{M_{X}}{\;\mathrm{TeV}}:\left\\{\begin{array}[]{l}\displaystyle\left[0.1:10\right]\quad({\rm
heavy}~{}M_{X})\\\\[8.53581pt]
\displaystyle\left[10^{-6}:0.08\right]\quad({\rm
light}~{}M_{X})\end{array}\right.\;.$ (43)

For the other three parameters, we use the same ranges for the two scenarios
of $M_{X}$ 555 There is also the possibility of both $M_{X}$ and $v_{\Phi}$
are light, which may lead to $Z_{3}=Z$ and both $Z_{1}$ and $Z_{2}$ are
lighter than $Z$. We will reserve this interesting scenario in future work.,

$\displaystyle 10^{-8}\leq$ $\displaystyle g_{H}$ $\displaystyle\leq
g^{\texttt{SM}}=\frac{e^{\texttt{SM}}}{\sin\theta_{W}}=0.65\;,$ $\displaystyle
10^{-8}\leq$ $\displaystyle g_{X}$ $\displaystyle\leq
g^{\prime\texttt{SM}}=\frac{e^{\texttt{SM}}}{\cos\theta_{W}}=0.35\;,$ (44)
$\displaystyle 5\;\mathrm{TeV}\leq$ $\displaystyle v_{\Phi}$
$\displaystyle\leq 200\;\mathrm{TeV}\;.$

We perform random scans by using MultiNest v2.17 Feroz:2008xx  with 30000
living points, an enlargement factor reduction parameter 0.5 and a stop
tolerance factor $10^{-3}$ . For sampling coverage, we combined several scans
and finally obtained $\sim 10^{5}$ samples for each scenario.

###  IV.2 Heavy $M_{X}$ Scenario

(a)

(b)

Figure 2:  Scatter plots in 1$\sigma$ on (a) ($M_{Z_{2}}$, ${\cal
O}^{2}_{22}$) plane and (b) ($M_{Z_{2}}$, ${\cal O}^{2}_{12}$) plane for the
heavy $M_{X}$ scenario. The red cross region with ${\cal O}_{22}^{2}$ between
0.8 and 1.0 represents the points of $W^{\prime 3}$-like $Z_{2}$ boson; the
blue triangle region with ${\cal O}_{22}^{2}$ between 0.2 and 0.8 represents
the points mixed with $W^{\prime 3}$, and the green circle region with ${\cal
O}_{22}^{2}$ between 0.0 and 0.2 represents the points of $X$-like $Z_{2}$
boson.

In the heavy $M_{X}$ scenario, the mass of $Z_{1}$ boson is located at around
$Z$-pole ($\sim 91\;\mathrm{GeV}$) so that $Z_{1}$ is identified as the SM
$Z$-boson. Note that $Z_{1}(Z)$ boson physics is strongly affected by the
different composition of $Z_{2}$ ($Z^{\prime}$) but not the heaviest boson
$Z_{3}$ ($Z^{\prime\prime}$) because $Z_{3}$ is heavier than $Z_{2}$ in our
parameter choices and therefore has less impact.

In Fig. 2, we present the scatter points of the composition of
$Z_{2}=\mathcal{O}_{12}Z^{\rm SM}+\mathcal{O}_{22}W^{\prime
3}+\mathcal{O}_{32}X$ for the $1\sigma$ region based on the likelihoods
described in Sec. IV.1. The color code hereafter represents the three
different composition of $Z_{2}$. Recalling Eq. (10), we define $W^{\prime
3}$-like $Z_{2}$ with condition $\mathcal{O}_{22}^{2}>0.8$ (red crosses
$\times$), mixed state $Z_{2}$ with $0.2<\mathcal{O}_{22}^{2}<0.8$ (blue
triangles $\triangledown$), and $X$-like $Z_{2}$ with condition
$\mathcal{O}_{22}^{2}<0.2$ (green circles $\circ$).

The $1\sigma$ allowed scatter points projected on the ($M_{Z_{2}}$,
$O^{2}_{22}$) and ($M_{Z_{2}}$, $O^{2}_{12}$) planes are depicted in Figs. 2a
and 2b, respectively. From the density of distribution in Fig. 2a, we can
clearly see that the mixed state $Z_{2}$ (blue triangles) is less evenly
distributed because it needs some trade-off between the two new gauge
couplings $g_{H}$ and $g_{X}$. In Fig. 2b, we projected the same parameter
space on the plane ($M_{Z_{2}}$, $\mathcal{O}_{12}^{2}$). Note that the mixing
$\mathcal{O}_{12}^{2}$ presents how $Z_{2}$ is consisted of $Z^{\rm{SM}}$.
Therefore, very small $\mathcal{O}_{12}^{2}$ implies
$\mathcal{O}_{32}^{2}\approx(1-\mathcal{O}_{22}^{2})$ from the orthogonality
of $\mathcal{O}$. Furthermore, the upper limit of $v_{\Phi}$ sets an lower
limit of the $\mathcal{O}_{12}^{2}$ for the red cross region. If $v_{\Phi}$
goes to infinity, $Z_{2}$ becomes super heavy and decouple. The composition of
$Z^{\rm{SM}}$ in $Z_{2}$ should then be negligible, thus
$\mathcal{O}_{12}^{2}$ vanishes in this limit. We note that the excluded
concave up region of $M_{Z_{2}}$ between 250 GeV and 6 TeV on the upper limit
of ${\cal O}_{12}^{2}$ is due to the constraint from ATLAS $Z^{\prime}$
search.

(a)

(b)

Figure 3: Scatter plots in $1\sigma$ on (a) ($M_{Z_{2}}$, $g_{H}$) plane and
(b) ($M_{Z_{2}}$, $v_{\Phi}$) plane for the heavy $M_{X}$ scenario. The color
code is the same as Fig. 2. The $1\sigma$ and $2\sigma$ contours of the
profile likelihood are also shown.

In Fig. 3, we show the $1\sigma$ (dashed) and $2\sigma$ (solid) likelihood
contours with scatter points inside the $1\sigma$ region on the (a)
($M_{Z_{2}}$, $g_{H}$) and (b) ($M_{Z_{2}}$, $v_{\Phi}$) planes. In Fig. 3a,
we can see that the $W^{\prime 3}$-like red crosses form a band with a
tendency proportional to $g_{H}$. This is because for a $W^{\prime 3}$-like
$Z_{2}$, $m^{2}_{Z_{2}}\approx g_{H}^{2}(v^{2}+v^{2}_{\Phi})/4\approx
g_{H}^{2}v^{2}_{\Phi}/4$ which can be extracted from the (3,3) element of the
mass matrix in Eq. (5). We can also see that at the lower bound of this band,
the 95% and 68% C.L. contours are overlapped because this lower bound is due
to our choice of $v_{\Phi}<200\;\mathrm{TeV}$ in its upper scan range, not
from the likelihood results. This implies that in the upper edge of this red
band where $g_{H}$ has larger value, the value of $v_{\Phi}$ there is smaller.
Therefore, the upper bound of this red cross band corresponds to the lower
values of $v_{\Phi}$, which can be excluded by the $\chi^{2}$ tolerance as we
can see in Fig. 3b where the scatter plot is projected on the ($M_{Z_{2}}$,
$v_{\Phi}$) plane. Surprisingly, in Fig. 3a, the blue triangle band,
corresponding to mixing mostly between $W^{\prime 3}$ and $X$ bosons, matches
the red cross band. This can be understood as the mass of $Z_{2}$ being
dominated by the (3,3) element of Eq. (5) even for an $80\%$ $X$ boson
composition. In the same figure, we can see the green circles running from
below the two red cross and blue triangle bands up to the upper limit of
$g_{H}$. In other words, we can see how the $M_{Z_{2}}$ passes from being
dominated by the (3,3) element of Eq. (5) (red crosses), which is
$g_{H}$-dependent, to being dominated by the $g_{H}$-independent (4,4) element
(green circles).

One particular feature of Fig. 3b is that the low $M_{Z_{2}}$ and low
$v_{\Phi}$ region (lower left corner) is covered only with $X$-like points
while both $W^{\prime 3}$-like and mixed points only approach this corner down
to a curved bound. This curved section in the lower bound can be related to
the curved upper bound for $W^{\prime 3}$ and mixed points in Fig. 3a for low
$M_{Z_{2}}<200$ GeV and $g_{H}\lesssim 10^{-2}$. These curves in the upper
bound (Fig. 3a) and in the lower bound (Fig. 3b) can be understood as smaller
$v_{\Phi}$ requiring larger $g_{H}$ to pass EWPT. In particular, if $g_{H}$ is
small, $v_{\Phi}$ has to be large in order to have a sizable diagonal (3,3)
element in the mass matrix in Eq. (5), while the off-diagonal (2,3) and (3,2)
elements remain small. However, the mixing effects from the off-diagonal
elements are not negligible and expected to be stronger when the ${Z_{2}}$
mass is getting closer to the ${Z^{\text{SM}}}$ mass. This gives rise to the
upper and lower bounds that we see in Figs. 3a and 3b, respectively, for the
$W^{\prime 3}$-like points. Such behaviour is not displayed for the $X$-like
points since they do not depend strongly on $g_{H}$.

The ATLAS $Z^{\prime}$ constraint almost rules out the region
$250\;\mathrm{GeV}<M_{Z_{2}}<6\;\mathrm{TeV}$ for $W^{\prime 3}$-like and
mixed $Z_{2}$, except the region with $v_{\Phi}>100\;\mathrm{TeV}$. However,
the $X$-like $Z_{2}$ at the same region has not been affected much by the
ATLAS $Z^{\prime}$ constraint.

Similarly, in Fig. 4, we show the $1\sigma$ (dashed) and $2\sigma$ (solid)
likelihood contours with scatter points inside the $1\sigma$ region on the (a)
($M_{Z_{3}}$, $g_{H}$) and (b) ($M_{Z_{3}}$, $g_{X}$) planes. From Fig. 4a,
one can easily see that the $X$-like $Z_{2}$ boson (green circles) forms a
band whose tendency is proportional to the $g_{H}$. This can be understood by
the fact that the composition of the $Z_{3}$ in this case is mainly from
$W^{\prime 3}$, which has a mass proportional to
$0.5\,g_{H}\sqrt{v^{2}+v_{\Phi}^{2}}\approx 0.5\,g_{H}v_{\Phi}$ again coming
from the (3,3) element of the mass matrix in Eq. (5). On the other hand, in
the case of the $W^{\prime 3}$-like $Z_{2}$ boson (red crosses), the mass of
the $Z_{3}$ almost does not depend on $g_{H}$. Indeed, the composition of
$Z_{3}$ is now mainly from $X$ and
$M^{2}_{Z_{3}}\approx(g_{X}^{2}(v^{2}+v_{\Phi}^{2})+M_{X}^{2})$. This is
clearly shown in Fig. 4b, when $g_{X}$ is small ($g_{X}<3\times 10^{-3}$), the
mass of $Z_{3}$ in the red cross region is dominated by $M_{X}$ and less than
our set-up limit of $10^{4}$ GeV. However, when $g_{X}$ is getting bigger, the
mass of the $Z_{3}$ can be dominated by the $g_{X}v_{\Phi}$ term for
sufficiently large value of $v_{\Phi}$. We can also see that the EWPT data
sets upper bounds on $g_{H}$ and $g_{X}$. The excluded concave up region of
$250\;\mathrm{GeV}<M_{Z_{2}}<6\;\mathrm{TeV}$ in Fig. 4a for the $W^{\prime
3}$-like and mixed composition of $Z_{2}$ is again due to the ATLAS
$Z^{\prime}$ search which does not apply for the $X$-like case. As a result,
the ATLAS $Z^{\prime}$ search cannot constrain on $g_{X}$ for $W^{\prime
3}$-like points as clearly shown in Fig. 4b.

(a)

(b)

Figure 4: Scatter plots in $1\sigma$ on (a) ($M_{Z_{3}}$, $g_{H}$) plane and
(b) ($M_{Z_{3}}$, $g_{X}$) plane for the heavy $M_{X}$ scenario. The color
code is the same as Fig. 2. The $1\sigma$ and $2\sigma$ contours of the
profile likelihood are also shown.

###  IV.3 Light $M_{X}$ Scenario

(a)

(b)

Figure 5:  Scatter plots in $1\sigma$ region on (a) ($M_{Z_{3}}$,
$O^{2}_{23}$) plane and (b) ($M_{Z_{3}}$, $O^{2}_{13}$) plane for the light
$M_{X}$ scenario. The red cross represents the points of $W^{\prime 3}$-like
$Z_{3}$ boson, the blue triangle represents the points mixed states ($Z^{\rm
SM}$, $W^{\prime 3}$ and $X$) $Z_{3}$ boson, and the green circles represents
$X$-like $Z_{3}$ boson.

In the light $M_{X}$ scenario, we require that the mass of $Z_{2}$ boson is
always at around $Z$-pole ($\sim 91\;\mathrm{GeV}$). In this scenario, the
lightest $Z_{1}$ with mass less than the $Z$-boson mass can be the dark photon
or dark $Z$, while the conventional $Z^{\prime}$ is the heaviest boson
$Z_{3}$. We note that the composition of $Z_{3}$ is given by
$Z_{3}=\mathcal{O}_{13}Z^{\rm SM}+\mathcal{O}_{23}W^{\prime
3}+\mathcal{O}_{33}X$. The $1\sigma$ allowed scatter points projected on the
($M_{Z_{3}}$, $O^{2}_{23}$) and ($M_{Z_{3}}$, $O^{2}_{13}$) planes are
depicted in Figs. 5a and 5b, respectively. The color code for the composition
of $Z_{3}$ is the same as in Fig. 2 for $Z_{2}$.

An obvious feature of Fig. 5a is that the mixed state of $Z_{3}$ (blue
triangles) has a mass upper limit. Intuitively, it requires some trade-off
between the gauge couplings $g_{H}$ and $g_{X}$ which results in
$M_{Z_{3}}\lesssim 500\;\mathrm{GeV}$. This effect will be discussed with more
detail later in Fig. 6. In Fig. 5b, we can see that the $Z^{\rm{SM}}$
composition of $Z_{3}$ is again small. However, unlike the heavy $M_{X}$
scenario, the $X$-like $Z_{3}$ boson has a similar distribution as $W^{\prime
3}$-like $Z_{3}$ boson. Additionally, the mixed $Z_{3}$ state at the mass
region between $210\;\mathrm{GeV}$ and $700\;\mathrm{GeV}$ cannot be excluded
by the ATLAS $Z^{\prime}$ constraint which is also different from the heavy
$M_{X}$ scenario.

(a)

(b)

Figure 6: Scatter plots in $1\sigma$ on (a) ($g_{H}$ $M_{Z_{3}}$) plane and
(b) ($v_{\Phi}$, $M_{Z_{3}}$) plane for the light $M_{X}$ scenario. The
markers are the same as Fig. 5. The $1\sigma$ and $2\sigma$ contours of the
profile likelihood are also shown.

In analogous to Fig. 3, we show in Fig. 6 the $1\sigma$ (dashed) and $2\sigma$
(solid) likelihood contours with scatter points in the $1\sigma$ region on the
(a) ($M_{Z_{3}}$, $g_{H}$) and (b) ($M_{Z_{3}}$, $v_{\Phi}$) planes. Comparing
Figs. 3a and 6a, we have a clear separation between the $W^{\prime 3}$-like
(red crosses) and $X$-like (green circles) regions in this light $M_{X}$
scenario. As before, the $W^{\prime 3}$-like red crosses follow a tendency
proportional to $g_{H}$ again because of the dominance of the (3,3) element of
Eq. (5) in $M_{Z_{3}}$, i.e., $M_{Z_{3}}\approx g_{H}v_{\Phi}/2$. Other
features shared between $W^{\prime 3}$-like points in Figs. 3a and 6a are the
distribution of $v_{\Phi}$ values; the $g_{H}$ lower bound owes to $v_{\Phi}$
upper bound but its upper bound owes to $\chi^{2}$ tolerance. As expected, the
ATLAS $Z^{\prime}$ search can constrain $g_{H}$ and $v_{\Phi}$ at the mass
region $250\;\mathrm{GeV}<M_{Z_{2}}<6\;\mathrm{TeV}$. However, the gauge
coupling for $X$-like $Z_{3}$ is proportional to $g_{X}$ not $g_{H}$ so that
the ATLAS $Z^{\prime}$ search cannot constrain on $g_{H}$ at the $X$-like
region, indicated by green circles. The $X$-like region in Fig. 6a has a
$g_{H}$ upper bound around $10^{-2}$ given by the $\chi^{2}$ tolerance and
most likely related to the lower bound on $v_{\Phi}$ displayed on Fig. 6b. The
mass of $Z_{3}$, $M_{Z_{3}}$, in this $X$-like green region can be
approximated by $\sqrt{g^{2}_{X}v^{2}_{\Phi}+M_{X}^{2}}$, this is why there is
not a clear $g_{H}$ dependence as in the $W^{\prime 3}$-like points. In Fig.
6a, as one would expect, the mixed region corresponds approximately to the
intersection between $X$-like and $W^{\prime 3}$-like regions, extending
lightly into their exclusive regions. This means that the upper and lower
bound of the mixed region are approximately given by the upper bound of the
$X$-like region and the lower bound of the $W^{\prime 3}$-region,
respectively. If we increase our maximum $v_{\Phi}$ value, the lower bound of
the $W^{\prime 3}$-like region would reach lower $g_{H}$ and the maximum
$M_{Z_{3}}$ for the mixed region would be increased. This is more clear after
looking at Fig. 6b where the maximum $M_{Z_{3}}$ value for the three regions
grows with the value of $v_{\Phi}$.

(a)

(b)

Figure 7: Scatter plots in $1\sigma$ on (a) ($g_{H}$ $M_{Z_{1}}$) plane and
(b) ($g_{X}$, $M_{Z_{1}}$) plane for the light $M_{X}$ scenario. The markers
are the same as Fig. 5. The $1\sigma$ and $2\sigma$ contours of the profile
likelihood are also shown.

Similarly, in Fig. 7, we show the $1\sigma$ (dashed) and $2\sigma$ (solid)
likelihood contours with scatter points in the $1\sigma$ region on the (a)
($M_{Z_{1}}$, $g_{H}$) plane and (b) ($M_{Z_{1}}$, $g_{X}$) planes. Again, the
red cross represents the points of $W^{\prime 3}$-like $Z_{3}$ boson, the blue
triangle represents the points of mixed state ($Z^{\rm SM}$, $W^{\prime 3}$
and $X$) $Z_{3}$ boson, and the green circle represents $X$-like $Z_{3}$
boson. We note that in this scenario, $Z_{1}$ is considered as a dark photon
666The lightest boson $Z_{1}$ can be tested in the dark photon experiments but
it is beyond the scope of this work. We will return to this in the future. and
has mass range from 1 MeV to $Z$-pole. One can easily see that the $Z_{3}$
composition is clearly separated on the planes of ($M_{Z_{1}},g_{H}$) and
($M_{Z_{1}},g_{X}$). In particular, while the $W^{\prime 3}$-like $Z_{3}$
boson parameter space is distributed in the region of larger $g_{H}$ and
smaller $g_{X}$, the $X$-like $Z_{3}$ boson, in contrast, prefers to be in the
region of smaller $g_{H}$ and larger $g_{X}$. The mixed composition of $Z_{3}$
lies in the range of $7\times 10^{-4}<g_{H}<5\times 10^{-3}$ and $4\times
10^{-4}<g_{X}<3\times 10^{-3}$. For the $X$-like $Z_{3}$ boson region in Fig.
7a, there is a lower bound for $g_{H}$ due to our choice of 200 TeV as the
upper bound for $v_{\Phi}$. Moreover, in Fig. 7b, one can also see that the
$\chi^{2}$ tolerance sets an upper limit on $g_{X}$ as the $Z_{1}$ boson mass
gets heavier.

Finally, we would like to emphasize that the contact interaction exclusion
regions at $M_{Z_{1}}<200\;\mathrm{MeV}$ and $10^{-4}<g_{X},g_{H}<10^{-3}$ are
owing to two different coupling components, $g_{H}\mathcal{O}_{2i}$ and
$g_{X}\mathcal{O}_{3i}$ in Eqs. (17) and (18).

###  IV.4 Future Prospects

Since current LEP together with other constraints already put a severe limit
on the parameter space, it will be interesting to see whether the future
$Z$-boson precision experiments can further probe our model. In the near
future, there are three colliders that can improve $Z$-boson measurements:
CEPC CEPC-SPPCStudyGroup:2015csa , ILC Fujii:2017vwa , and FCC-ee
dEnterria:2016fpc . Among them, CEPC is the one that could give the most
sensitive limit. Therefore, in this subsection, we make an estimation of our
parameter space with the projected CEPC sensitivity.

In the third column of Table 2, we quote the expected CEPC sensitivity CEPC-
SPPCStudyGroup:2015csa . Apparently, some of the error bars are expected to be
significantly reduced. Note that the CEPC preliminary conceptual design report
does not provide a full list as the LEP measurements showed in the 2nd column.
Therefore, for those missing rows, we reuse the data from the 2nd column (LEP
data).

Figure 8: The $\Delta\chi^{2}$ as function of $v_{\Phi}$. The red solid line
and blue dashed line are based on present constraint and future CEPC
sensitivity. The left and right panels are corresponding to heavy and light
$M_{X}$ scenarios, respectively.

To start with, we present the $\Delta\chi^{2}$ in terms of $v_{\Phi}$ in Fig.
8 for heavy (left) and light (right) $M_{X}$ scenarios. Importantly,
$v_{\Phi}$ is the most sensitive parameter in the G2HDM, determining the
theory scale. For the heavy $M_{X}$ scenario, in the present sensitivity case
the $2\sigma$ lower bound is around $24\;\mathrm{TeV}$, while in the CEPC case
it is around 44 TeV. For the light $M_{X}$ scenario, the 2$\sigma$ current and
CEPC lower limit of $v_{\Phi}$ is smaller than the heavier $M_{X}$ scenario.
In particular, $v_{\Phi}>15\text{~{}TeV}\,(36\text{~{}TeV})$ at $95\%$ C.L.
from current experiments (CEPC). The difference between these two scenarios is
owing to the different sources of constraints on $v_{\Phi}$. For the heavy
$M_{X}$ scenario, the EWPT constraints of the SM $Z$ boson play an important
role in raising the lower limit of $v_{\Phi}$. However, for the light $M_{X}$
scenario, the main constraint to exclude the lower $v_{\Phi}$ region is from
$Z^{\prime}$ searches. This also explains why the future sensitivity does not
further push $v_{\Phi}$ in the light $M_{X}$ scenario to larger values as the
heavy $M_{X}$ scenario does because the future sensitivities of contact
interactions are not available for CEPC and only the previous limits from LEP
II are used.

In Fig. 9, we compare the present limit and future CEPC sensitivity of the
two-dimensional contours on the ($g_{H},~{}g_{X}$) plane. The figure in the
left (right) column corresponds to the heavy (light) $M_{X}$ scenario. Because
the upper scan limit of $g_{H}$ is set to be less than $g_{2}^{\texttt{SM}}$,
the experimental constraints on $g_{H}$ are not present. In contrast, $g_{X}$
has an upper limit from the constraints due to $g_{H}$ having a lower limit
from the maximum scanned $v_{\Phi}$ value. The upper limit of $g_{H}$ can be
further improved by future CEPC sensitivity along the edge of the contour.
However, the light $M_{X}$ scenario is mildly constrained by future CEPC
sensitivity. The two contour plots in Fig. 9 can be further understood as
follows. We note that, for the case of $W^{\prime 3}$-like $Z_{2}$ in heavy
$M_{X}$ scenario (left panel) or $W^{\prime 3}$-like $Z_{3}$ in light $M_{X}$
scenario (right panel), $g_{H}$ has a lower limit at $\sim 2\times 10^{-3}$
due to our choices of the parameter scan ranges. Indeed, in both cases we have
$M_{Z_{2,3}}\approx 0.5g_{H}\sqrt{v^{2}+v_{\Phi}^{2}}$ which implies that
$g_{H}\approx(2M_{Z_{2,3}})/\sqrt{v^{2}+v_{\Phi}^{2}}$. Since we require
$M_{Z_{2,3}}>210$ GeV and $v_{\Phi}<200$ TeV, this implies $g_{H}>2\times
10^{-3}$. Similarly, for the case of $X$-like $Z_{3}$ in the light $M_{X}$
scenario (upper right panel), the mass of $Z_{3}$, is given by
$M_{Z_{3}}\approx\sqrt{g_{X}^{2}(v^{2}+v_{\Phi}^{2})+M_{X}^{2}}$ so that we
can obtain
$g_{X}\approx\sqrt{M_{Z_{3}}^{2}-M_{X}^{2}}/\sqrt{v^{2}+v_{\Phi}^{2}}$. This
yields a lower limit for $g_{X}$ at $\sim 10^{-3}$ when we require
$M_{Z_{3}}>210$ GeV, $M_{X}<80$ GeV and $v_{\Phi}<200$ TeV. On the other hand,
$g_{X}$ has no lower limit in the heavy $M_{X}$ scenario (left panel).

Figure 9: The present and future sensitivity allowed regions projected on the
($g_{H},~{}g_{X}$) plane in both heavy (left) and light (right) $M_{X}$
scenarios. The red solid line is for the present $95\%$ limit while the blue
dashed line is for the CEPC future sensitivity. Scatter points in $1\sigma$
region of the present constraint are also shown. The color codes in the left
and right panel are same as in Figs. 2 and Fig. 5, respectively.

The Stüeckelberg mass parameter $M_{X}$ is a filter to split the parameter
space into two scenarios but we have not been able to constrain this
parameter. The reason is simply that $Z_{3}$ in the heavy $M_{X}$ scenario is
too heavy to be relevant by current experiments. On the other hand, in the
light $M_{X}$ scenario with $M_{X}<80\;\mathrm{GeV}$, it is again too light to
be presented in the EWPT data. To constrain light $M_{X}$, just like dark
photon, the lightest $Z_{1}$ could be detected by those future beam dump
experiments such as NA62 1510.00172 , Belle II 1002.5012 , and SHiP 1504.04855
. However, this is beyond the scope of this work and we will return to it in
the future.

##  V Summary and Conclusion

In this paper, we perform an updated profile likelihood analysis for the gauge
sector of G2HDM.

For the two Stüeckelberg mass parameters $M_{Y}$ and $M_{X}$ associated with
the hypercharge $U(1)_{Y}$ and the extra $U(1)_{X}$ respectively, we showed
that a nonzero $M_{Y}$ would produce non-standard QED couplings for all the
fermions in G2HDM, albeit we can always achieve a massless photon for
arbitrary values of $M_{X}$ and $M_{Y}$. We therefore set $M_{Y}=0$ in our
numerical analysis. The remaining new parameters in the gauge sector of G2HDM
needed to be constrained are $g_{H}$, $g_{X}$, $v_{\Phi}$ and $M_{X}$.

We have examined the remaining parameter space with the EWPT LEP data at the
$Z$-pole, contact interaction constraints from LEP-II and LHC Run II data for
the search of high-mass dilepton resonances. The contact interactions
constraints can definitely provide a lower limit on $v_{\Phi}$, but the EWPT
data play a significant role to constrain the parameter space non-trivially.
While the LHC search for the high-mass dilepton resonances also impose
important constraints on the parameter space, the Drell-Yan data from the $Z$
decay does not impose noticeable impacts yet.

We classify our parameter space based on three different composition
($X$-like, $W^{\prime 3}$-like, and mixed) of the heavy neutral gauge boson,
either $Z_{2}$ or $Z_{3}$, which is the next-heavier $Z$ boson than the SM
one, in order to manifest the physics and constraints discussed in this paper.

In the heavy $M_{X}$ scenario ($M_{X}>100\;\mathrm{GeV}$), the SM-like $Z$ is
the lightest $Z_{1}$ boson and EWPT constraints exclude the small $v_{\Phi}$
region up to $24\;\mathrm{TeV}$ at $2\sigma$ significance. However, the EWPT
constraints are not so sensitive to the light $M_{X}$ scenario
($M_{X}<80\;\mathrm{GeV}$) where SM-like $Z$ is the next-lightest $Z_{2}$
boson. In particular, the $v_{\Phi}$ is required to be greater than
$15\;\mathrm{TeV}$ due to the constraints of $Z^{\prime}$ contact interaction
search from LEP-II and high-mass dilepton resonance search from LHC Run II.
Furthermore, in both light and heavy $M_{X}$ scenarios, $M_{X}$ is just a
parameter to tweak between two scenarios and it is totally unbounded in this
study. It is likely that the future dark photon searches might set a limit on
the $M_{X}$ in the light $M_{X}$ scenario. On the other hand, it is not so
trivial for the couplings $g_{X}$ and $g_{H}$ because we found it is hard to
set an upper bound on them individually.

Although the SM $Z$ boson is fixed at the $Z$-pole, the allowed physical
masses of the heavier $Z_{i}$ still depend on the $M_{X}$ and detailed
composition. Generally speaking, the $Z_{2}$ allowed mass range in the heavy
$M_{X}$ scenario is same as the range of $M_{X}$ but $Z_{3}$ mass can reach up
to $70\;\mathrm{TeV}$ for $X$-like composition and $40\;\mathrm{TeV}$ for both
$W^{\prime 3}$-like and mixed composition. Like the role of $M_{X}$ in the
heavy scenario, the $M_{Z_{1}}$ in the light $M_{X}$ scenario is dominated by
$M_{X}$ and the allowed mass ranges of $M_{Z_{1}}$ have no difference between
different composition. However, regarding to $M_{Z_{3}}$, mixed $Z_{3}$ is
restricted to less than $500\;\mathrm{GeV}$ but the masses of $X$-like and
$W^{\prime 3}$-like $Z_{3}$ are below $70\;\mathrm{TeV}$.

Finally, we also discuss the future sensitivity of the new parameters at the
CEPC. We found that the CEPC can significantly probe the parameter space of
the heavy $M_{X}$ scenario but the sensitivity is not improved much for the
light $M_{X}$ scenario. In the latter case, when $M_{X}$ is getting very
light, $Z_{1}$ can be much lighter than the $Z$-boson and it is more
appropriate to identify it as the dark photon or dark $Z_{D}$. The very rich
phenomenology of light dark photon or dark $Z_{D}$ in G2HDM remains to be
explored in the future.

## Acknowledgments

We would like to thank Wei-Chih Huang, Zuowei Liu and Xun Xue for stimulating
discussions. This work was supported in part by the Ministry of Science and
Technology (MoST) of Taiwan under Grant Nos. 107-2119-M-001-033- and
107-2811-M-001-027-. Y.-L. S. Tsai was funded in part by Chinese Academy of
Sciences Taiwan Young Talent Programme under Grant No. 2018TW2JA0005. V. Q.
Tran was funded in part by the National Natural Science Foundation of China
under Grant Nos. 11775109, U1738134 and by the National Recruitment Program
for Young Professionals.

##  Appendix A The Rotation Angles $\phi$, $\theta$ and $\psi$

In this appendix, we will show how to obtain the equations of the rotation
angles such as Eqs. (7), (8) and (9) from the orthogonal matrix which
diagonalizes the mass matrix $\mathcal{M}^{2}_{\rm gauge}(M_{Y}=0)$ given in
Eq. (5). The orthogonal matrix we choose is Eq. (6) because it is rather
convenient to find all the ${\cal O}_{ij}$s and determine the rotational
angles $\phi$, $\theta$ and $\psi$ numerically. However, the computation of
the angles in terms of the fundamental parameters in the Lagrangian are
difficult to organize into nice forms using Eq. (6) for $\mathcal{O}$, so we
apply Cramer’s rule for solving the secular equations and get another form for
$\mathcal{O}$ as follows

$\mathcal{O}\ =\
\left(\begin{array}[]{ccc}|x_{1}|/\Delta_{1}&x_{2}/\Delta_{2}&x_{3}/\Delta_{3}\\\
y_{1}/\Delta_{1}&|y_{2}|/\Delta_{2}&y_{3}/\Delta_{3}\\\
z_{1}/\Delta_{1}&z_{2}/\Delta_{2}&|z_{3}|/\Delta_{3}\end{array}\right)\,,$
(45)

where

$\Delta_{i}=\sqrt{x^{2}_{i}\,+\,y^{2}_{i}\,+\,z^{2}_{i}}\,,$ (46)

and

$\displaystyle
x_{1}=\\!\left|\begin{array}[]{lr}M^{2}_{22}-M^{2}_{Z_{1}}&M^{2}_{23}\\\
M^{2}_{32}&M^{2}_{33}-M^{2}_{Z_{1}}\end{array}\right|\\!\,,\ y_{1}={\rm
s}_{x_{1}}\\!\left|\begin{array}[]{cc}M^{2}_{23}&M^{2}_{21}\\\
M^{2}_{33}-M^{2}_{Z_{1}}&M^{2}_{31}\end{array}\right|\,,\ z_{1}={\rm
s}_{x_{1}}\\!\left|\begin{array}[]{cc}M^{2}_{21}&M^{2}_{22}-M^{2}_{Z_{1}}\\\
M^{2}_{31}&M^{2}_{32}\end{array}\right|\,,$ (53) $\displaystyle x_{2}={\rm
s}_{y_{2}}\\!\left|\begin{array}[]{cc}M^{2}_{13}&M^{2}_{12}\\\
M^{2}_{33}-M^{2}_{Z_{2}}&M^{2}_{32}\end{array}\right|\,,\
y_{2}=\\!\left|\begin{array}[]{lr}M^{2}_{11}-M^{2}_{Z_{2}}&M^{2}_{13}\\\
M^{2}_{31}&M^{2}_{33}-M^{2}_{Z_{2}}\end{array}\right|\\!\,,\ z_{2}={\rm
s}_{y_{2}}\\!\left|\begin{array}[]{cc}M^{2}_{12}&M^{2}_{11}-M^{2}_{Z_{2}}\\\
M^{2}_{32}&M^{2}_{31}\end{array}\right|\,,$ (60) $\displaystyle x_{3}={\rm
s}_{z_{3}}\\!\left|\begin{array}[]{cc}M^{2}_{12}&M^{2}_{13}\\\
M^{2}_{22}-M^{2}_{Z_{3}}&M^{2}_{23}\end{array}\right|,\ y_{3}={\rm
s}_{z_{3}}\\!\left|\begin{array}[]{cc}M^{2}_{13}&M^{2}_{11}-M^{2}_{Z_{3}}\\\
M^{2}_{23}&M^{2}_{21}\end{array}\right|\,,\
z_{3}=\\!\left|\begin{array}[]{lr}M^{2}_{11}-M^{2}_{Z_{3}}&M^{2}_{12}\\\
M^{2}_{21}&M^{2}_{22}-M^{2}_{Z_{3}}\end{array}\right|\\!\,,$ (67)

with $M^{2}_{ij}$ stands for the element of $\mathcal{M}^{2}_{\rm
gauge}(M_{Y}=0)$, $M^{2}_{Z_{i}}(i=1,2,3)$ are the mass eigenvalues and ${\rm
s}_{x_{i}}={\rm sign}(x_{i})$. From Eq. (6), one can obtain the following
relations for the rotational angles $\phi$, $\theta$ and $\psi$ 777We note
that similar approach had been used in Pilaftsis:1999qt  for the scalar boson
mass matrix in MSSM with explicit CP violation.

$\phi=\arctan{\left(\frac{-\mathcal{O}_{12}}{\mathcal{O}_{22}}\right)},\,\theta=\arctan\left({\frac{-\mathcal{O}_{32}}{\mathcal{O}_{12}}}\sin{\phi}\right),\,\psi={\rm
arccot}\left({\frac{-\mathcal{O}_{21}}{\mathcal{O}_{31}}}\frac{\cos{\theta}}{\sin{\phi}}-\sin{\theta}\cot{\phi}\right)\,,$
(69)

with the range for $\theta$ covers $\pi$ radians, and the range for $\phi$ and
$\psi$ covers $2\pi$ radians. Note that the expressions in Eq. (69) do not
depend on the $\Delta_{i}$ given in Eq. (46). Using Eqs. (45) and (53) for the
various $\mathcal{O}_{ij}$ in Eq. (69), after some algebra, one can obtain
Eqs. (7), (8) and (9), which are collected here again for convenience.

$\tan({\phi})=\frac{-g_{H}vM_{Z^{\text{SM}}}(M_{X}^{2}-M_{Z_{2}}^{2}+2g_{X}^{2}v_{\Phi}^{2})}{2\Big{(}M_{Z_{2}}^{4}-\big{(}M_{Z^{\text{SM}}}^{2}+M_{X}^{2}+(v^{2}+v_{\Phi}^{2})g_{X}^{2}\big{)}M_{Z_{2}}^{2}+M_{Z^{\text{SM}}}^{2}(M_{X}^{2}+g_{X}^{2}v_{\Phi}^{2})\Big{)}}\,,$
(70)
$\tan({\theta})=\frac{-g_{X}(M_{Z_{2}}^{2}(v^{2}-v_{\Phi}^{2})+M_{Z^{\text{SM}}}^{2}v_{\Phi}^{2})}{vM_{Z^{\text{SM}}}(M_{X}^{2}-M_{Z_{2}}^{2}+2g_{X}^{2}v_{\Phi}^{2})}\sin{\phi}\,,$
(71)
$\cot({\psi})=\frac{g_{H}(M_{Z_{1}}^{2}-M_{X}^{2}-2g_{X}^{2}v_{\Phi}^{2})}{g_{X}(g_{H}^{2}v_{\Phi}^{2}-2M_{Z_{1}}^{2})}\frac{\cos{\theta}}{\sin{\phi}}-\sin{\theta}\cot{\phi}\,.$
(72)

Thus one can compute the rotation angles in terms of the fundamental
parameters of the model which can provide some useful insights in the
vanishing limits of $g_{H}$ and $g_{X}$ as discussed in Section (II.2).

##  Appendix B Decay Widths of New Neutral Gauge Bosons

(a)

(b)

Figure 10: Heavy $M_{X}$ scenario: scatter plots in $1\sigma$ region on the
(a) ($\Gamma_{Z_{2}}$/$M_{Z_{2}}$, $M_{Z_{2}}$) and (b)
($\Gamma_{Z_{3}}$/$M_{Z_{3}}$, $M_{Z_{3}}$) planes. The markers are the same
as Fig. 2.

(a)

(b)

Figure 11: Light $M_{X}$ scenario: scatter plots in $1\sigma$ on the (a)
($\Gamma_{Z_{3}}$/$M_{Z_{3}}$, $M_{Z_{3}}$) and (b)
($\Gamma_{Z_{1}}$/$M_{Z_{1}}$, $M_{Z_{1}}$) planes. The markers are the same
as Fig. 5.

In this appendix, we show the decay widths of the two new neutral gauge bosons
$Z_{i}$. We note that for light $M_{X}$ scenario, $i={(1,3)}$, while for heavy
$M_{X}$ scenario, $i={(2,3)}$.

  * •

The decay width of $Z_{i}$ to a pair of fermions (including both SM and new
heavy fermions) is given as follows

$\Gamma(Z_{i}\rightarrow
f\bar{f})=\frac{N^{c}_{f}g_{M}^{2}M_{Z_{i}}}{12\pi}\sqrt{1-4r_{if}}\left((2r_{if}+1)|{v^{(i)}_{f}}|^{2}+(1-4r_{if})|{a^{(i)}_{f}}|^{2}\right)\,,$
(73)

where $g_{M}=\sqrt{g^{2}+g^{\prime 2}}/2$, $N^{c}_{f}$ is the number of color
for fermion $f$, the coefficients $v_{f}^{(i)}$ and $a_{f}^{(i)}$ are the
couplings that appear in Eqs. (17) and (18) and
$r_{if}=\frac{m_{f}^{2}}{M_{Z_{i}}^{2}}.$

  * •

The decay width for $Z_{i}\rightarrow W^{+}W^{-}$ process is given by
Barger:1987xw

$\Gamma(Z_{i}\rightarrow W^{+}W^{-})=\frac{g_{Z_{i}WW}^{2}M_{Z_{i}}}{192\pi
r_{iW}^{2}}(1-4r_{iW})^{3/2}\Big{(}1+20r_{iW}+12r_{iW}^{2}\Big{)}\,,$ (74)

where $r_{iW}=\frac{M_{W}^{2}}{M_{Z_{i}}^{2}}$ and the coupling
$g_{Z_{i}WW}=gc_{W}{\cal O}_{1i}$ .

  * •

Similarly, one can obtain the decay width for $Z_{i}\rightarrow W^{\prime
p}W^{\prime m}$ process as

$\Gamma(Z_{i}\rightarrow W^{\prime p}W^{\prime
m})=\frac{g_{Z_{i}W^{\prime}W^{\prime}}^{2}M_{Z_{i}}}{192\pi
r_{iW^{\prime}}^{2}}(1-4r_{iW^{\prime}})^{3/2}\Big{(}1+20r_{iW^{\prime}}+12r_{iW^{\prime}}^{2}\Big{)}\,,$
(75)

where $r_{iW^{\prime}}=\frac{M_{W^{\prime}}^{2}}{M_{Z_{i}}^{2}}$ and the
coupling $g_{Z_{i}W^{\prime}W^{\prime}}=g_{H}\,{\cal O}_{2i}$ .

  * •

The new neutral gauge boson $Z_{i}$ can also decay into pair of scalar dark
matter candidate in this model. The decay width for this process
$Z_{i}\rightarrow DD^{*}$ is given by Barger:1987xw

$\Gamma(Z_{i}\rightarrow
DD^{*})=\frac{g_{Z_{i}DD}^{2}M_{Z_{i}}}{48\pi}\,(1-4r_{iD})^{3/2}\,,$ (76)

where the coupling $g_{Z_{i}DD}=g_{H}{\cal O}_{2i}$ and
$r_{iD}=\frac{M_{D}^{2}}{M_{Z_{i}}^{2}}$ . We note that $D$ is a triplet-like
scalar dark matter in this model and we assumed this dark matter doesn’t mix
with other scalars in this calculation.

  * •

The decay width for $Z_{i}\rightarrow H^{+}H^{-}$ is given by

$\Gamma(Z_{i}\rightarrow
H^{+}H^{-})=\frac{g_{Z_{i}H^{+}H^{-}}^{2}M_{Z_{i}}}{48\pi}\,(1-4r_{iH^{\pm}})^{3/2}\,,$
(77)

where $r_{iH^{\pm}}=\frac{M_{H^{\pm}}^{2}}{M_{Z_{i}}^{2}}$ and the coupling
$g_{Z_{i}H^{+}H^{-}}$ is given as follows

$g_{Z_{i}H^{+}H^{-}}=\frac{1}{2}(c_{W}g-s_{W}g^{\prime}){\cal
O}_{1i}-\frac{1}{2}g_{H}{\cal O}_{2i}+g_{X}{\cal O}_{3i}\,.$ (78)

  * •

The decay width for $Z_{i}\rightarrow Z_{j}H$ is given by Barger:1987xw

$\displaystyle\Gamma(Z_{i}\rightarrow
Z_{j}H)=\frac{g_{Z_{i}Z_{j}H}^{2}M_{Z_{i}}}{192\pi\,M_{Z_{j}}^{2}}\,\Big{(}1+(r_{ij}-r_{iH})^{2}-2(r_{ij}+r_{iH})\Big{)}^{1/2}$
$\displaystyle\times\Big{(}1+(r_{ij}-r_{iH})^{2}+10r_{ij}-2r_{iH}\Big{)}\,,$
(79)

where $r_{iH}=\frac{M_{H}^{2}}{M_{Z_{i}}^{2}}$,
$r_{ij}=\frac{M_{Z_{j}}^{2}}{M_{Z_{i}}^{2}}$ and the coupling
$g_{Z_{i}Z_{j}H}$ is given as follows

$\displaystyle g_{Z_{i}Z_{j}H}=\frac{v}{2}\Big{(}(c_{W}g+s_{W}g^{\prime}){\cal
O}_{1j}-g_{H}{\cal O}_{2j}-2g_{X}{\cal O}_{3j}\Big{)}$
$\displaystyle\times\Big{(}(c_{W}g+s_{W}g^{\prime}){\cal O}_{1i}-g_{H}{\cal
O}_{2i}-2g_{X}{\cal O}_{3i}\Big{)}\,,$ (80)

here $v$ is the VEV of the SM Higgs field. Note that we have ignored the
mixing of SM like-Higgs $H$ with other scalar bosons in the above
calculations.

  * •

Finally, if not kinematically prohibited, the new neutral gauge bosons can
also decay into $W^{\prime}$ and the dark matter $D$. The decay width for this
process can be computed as

$\displaystyle\Gamma(Z_{i}\rightarrow W^{\prime p}D^{*}/W^{\prime
m}D)=\frac{g_{Z_{i}W^{\prime}D}^{2}M_{Z_{i}}}{192\pi\,M^{2}_{W^{\prime}}}\,\Big{(}1+(r_{iW^{\prime}}-r_{iD})^{2}-2(r_{iW^{\prime}}+r_{iD})\Big{)}^{1/2}$
$\displaystyle\times\Big{(}1+(r_{iW^{\prime}}-r_{iD})^{2}+10r_{iW^{\prime}}-2r_{iD}\Big{)}\,,$

where the coupling $g_{Z_{i}W^{\prime}D}=g_{H}^{2}{\cal O}_{2i}v_{\Delta}$
with $v_{\Delta}$ being the VEV of $SU(2)_{H}$ triplet Higgs.

In Figs. 10 and 11, we show the scatter plots of the ratio of decay width over
mass of the two new gauge bosons in the heavy $M_{X}$ and light $M_{X}$
scenarios respectively. In those plots, we set the dark matter mass $M_{D}$ to
be $10\%$ of the new heavy neutral gauge boson $Z_{i}$ ($i.e.$
$M_{D}=0.1\times M_{Z_{2}}$ in the case of heavy $M_{X}$ scenario, while
$M_{D}=0.1\times M_{Z_{3}}$ in the case of light $M_{X}$ scenario), charged
Higgs mass $M_{H^{\pm}}$ equals 1.5 TeV and the mass of $W^{\prime(p,m)}$ is
randomly chosen in the range of [$M_{D},200$ TeV]. Moreover, we assume the
masses of new heavy fermions are degenerate and equal to 3 TeV. We note that
$v_{\Delta}$ can be derived from other parameters according to
$v_{\Delta}=0.5\sqrt{\frac{4M_{W^{\prime}}^{2}}{g_{H}^{2}}-(v^{2}+v_{\Phi}^{2})}$.
From these scatter plots one can see that for the heavy neutral gauge bosons
in both scenarios, their ratios $\Gamma_{Z_{i}}/M_{Z_{i}}$ are all below $\sim
1\%$, until they are heavier than 10 TeV the ratios can then reach $\sim 6\%$.
However for the light $Z_{1}$ in the light $M_{X}$ scenario,
$\Gamma_{Z_{1}}/M_{Z_{1}}$ is well below $10^{-4}$.

## References

  * (1) A. Salam and J. C. Ward, “Weak and electromagnetic interactions,” Nuovo Cim. 11, 568 (1959). doi:10.1007/BF02726525
  * (2) S. L. Glashow, “Partial Symmetries of Weak Interactions,” Nucl. Phys. 22, 579 (1961). doi:10.1016/0029-5582(61)90469-2
  * (3) S. Weinberg, “A Model of Leptons,” Phys. Rev. Lett. 19, 1264 (1967). doi:10.1103/PhysRevLett.19.1264
  * (4) A. Salam (1968). N. Svartholm, ed. Elementary Particle Physics: Relativistic Groups and Analyticity. Eighth Nobel Symposium. Stockholm: Almquvist and Wiksell. p. 367.
  * (5) F. Englert and R. Brout, “Broken Symmetry and the Mass of Gauge Vector Mesons,” Phys. Rev. Lett. 13, 321 (1964). doi:10.1103/PhysRevLett.13.321
  * (6) P. W. Higgs, “Broken Symmetries and the Masses of Gauge Bosons,” Phys. Rev. Lett. 13, 508 (1964). doi:10.1103/PhysRevLett.13.508
  * (7) G. S. Guralnik, C. R. Hagen and T. W. B. Kibble, “Global Conservation Laws and Massless Particles,” Phys. Rev. Lett. 13, 585 (1964). doi:10.1103/PhysRevLett.13.585
  * (8) For a recent review on dark matter, see T. Lin, “TASI lectures on dark matter models and direct detection,” arXiv:1904.07915 [hep-ph] and references therein.
  * (9) For a review on dark energy, see for example P. J. E. Peebles and B. Ratra, “The Cosmological constant and dark energy,” Rev. Mod. Phys. 75, 559 (2003) doi:10.1103/RevModPhys.75.559 [astro-ph/0207347] and references therein.
  * (10) See for example the textbook by S. Dodelson, “Modern Cosmology,” Amsterdam, Netherlands: Academic Pr. (2003).
  * (11) W. C. Huang, Y. L. S. Tsai and T. C. Yuan, “G2HDM : Gauged Two Higgs Doublet Model,” JHEP 1604, 019 (2016) doi:10.1007/JHEP04(2016)019 [arXiv:1512.00229 [hep-ph]].
  * (12) N. G. Deshpande and E. Ma, “Pattern of symmetry breaking with two Higgs doublets,” Phys. Rev. D 18, 2574 (1978). doi:10.1103/PhysRevD.18.2574
  * (13) E. Ma, “Verifiable radiative seesaw mechanism of neutrino mass and dark matter,” Phys. Rev. D 73, 077301 (2006) doi:10.1103/PhysRevD.73.077301 [hep-ph/0601225].
  * (14) 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) doi:10.1103/PhysRevD.74.015007 [hep-ph/0603188].
  * (15) 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) doi:10.1088/1475-7516/2007/02/028 [hep-ph/0612275].
  * (16) 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. Rep. 516, 1 (2012) doi:10.1016/j.physrep.2012.02.002 [arXiv:1106.0034 [hep-ph]].
  * (17) 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 (2014) 030 doi:10.1088/1475-7516/2014/06/030 [arXiv:1310.0358 [hep-ph]].
  * (18) A. Arhrib, R. Benbrik and T. C. Yuan, “Associated Production of Higgs at Linear Collider in the Inert Higgs Doublet Model,” Eur. Phys. J. C 74, 2892 (2014) doi:10.1140/epjc/s10052-014-2892-5 [arXiv:1401.6698 [hep-ph]].
  * (19) A. Ilnicka, M. Krawczyk, and T. Robens, “Inert doublet model in light of LHC Run I and astrophysical data,” Phys. Rev. D 93, no. 5, 055026 (2016) doi:10.1103/PhysRevD.93.055026 [arXiv:1508.01671 [hep-ph]].
  * (20) 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) doi:10.1103/PhysRevD.97.035011 [arXiv:1612.00511 [hep-ph]].
  * (21) B. Eiteneuer, A. Goudelis, and J. Heisig, “The inert doublet model in the light of Fermi-LAT gamma-ray data: a global fit analysis,” Eur. Phys. J. C 77, no. 9, 624 (2017) doi:10.1140/epjc/s10052-017-5166-1 [arXiv:1705.01458 [hep-ph]].
  * (22) D. Borah, P. S. B. Dev and A. Kumar, “TeV scale leptogenesis, inflaton dark matter and neutrino mass in a scotogenic model,” Phys. Rev. D 99, no. 5, 055012 (2019) doi:10.1103/PhysRevD.99.055012 [arXiv:1810.03645 [hep-ph]].
  * (23) T. W. Kephart and T. C. Yuan, “Origins of inert Higgs doublets,” Nucl. Phys. B906, 549 (2016) doi:10.1016/j.nuclphysb.2016.03.023 [arXiv:1508.00673 [hep-ph]].
  * (24) A. Goudelis, B. Herrmann and O. Stål, “Dark matter in the Inert Doublet Model after the discovery of a Higgs-like boson at the LHC,” JHEP 1309, 106 (2013) doi:10.1007/JHEP09(2013)106 [arXiv:1303.3010 [hep-ph]].
  * (25) B. Swiezewska and M. Krawczyk, “Diphoton rate in the inert doublet model with a 125 GeV Higgs boson,” Phys. Rev. D 88, no. 3, 035019 (2013) doi:10.1103/PhysRevD.88.035019 [arXiv:1212.4100 [hep-ph]].
  * (26) A. Arhrib, R. Benbrik and N. Gaur, “$H\rightarrow\gamma\gamma$ in Inert Higgs Doublet Model,” Phys. Rev. D 85, 095021 (2012) doi:10.1103/PhysRevD.85.095021 [arXiv:1201.2644 [hep-ph]].
  * (27) L. M. Krauss and F. Wilczek, “Discrete Gauge Symmetry in Continuum Theories,” Phys. Rev. Lett. 62, 1221 (1989). doi:10.1103/PhysRevLett.62.1221
  * (28) R. Kallosh, A. D. Linde, D. A. Linde, and L. Susskind, “Gravity and global symmetries,” Phys. Rev. D 52, 912 (1995) doi:10.1103/PhysRevD.52.912 [hep-th/9502069].
  * (29) C. R. Chen, Y. X. Lin, C. S. Nugroho, R. Ramos, Y. L. S. Tsai and T. C. Yuan, in preparation.
  * (30) C. R. Chen, Y. X. Lin, V. Q. Tran and T. C. Yuan, “Pair production of Higgs bosons at the LHC in gauged 2HDM,” Phys. Rev. D 99, no. 7, 075027 (2019) doi:10.1103/PhysRevD.99.075027 [arXiv:1810.04837 [hep-ph]].
  * (31) 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) doi:10.1103/PhysRevD.98.095006 [arXiv:1806.05632 [hep-ph]].
  * (32) P. Ko, Y. Omura and C. Yu, “A Resolution of the Flavor Problem of Two Higgs Doublet Models with an Extra $U(1)_{H}$ Symmetry for Higgs Flavor,” Phys. Lett. B 717, 202 (2012) doi:10.1016/j.physletb.2012.09.019 [arXiv:1204.4588 [hep-ph]].
  * (33) M. D. Campos, D. Cogollo, M. Lindner, T. Melo, F. S. Queiroz and W. Rodejohann, “Neutrino Masses and Absence of Flavor Changing Interactions in the 2HDM from Gauge Principles,” JHEP 1708, 092 (2017) doi:10.1007/JHEP08(2017)092 [arXiv:1705.05388 [hep-ph]].
  * (34) D. A. Camargo, L. Delle Rose, S. Moretti and F. S. Queiroz, “Collider bounds on 2-Higgs doublet models with U (1)X gauge symmetries,” Phys. Lett. B 793, 150 (2019) doi:10.1016/j.physletb.2019.04.048 [arXiv:1805.08231 [hep-ph]].
  * (35) D. A. Camargo, A. G. Dias, T. B. de Melo and F. S. Queiroz, “Neutrino Masses in a Two Higgs Doublet Model with a U(1) Gauge Symmetry,” JHEP 1904, 129 (2019) doi:10.1007/JHEP04(2019)129 [arXiv:1811.05488 [hep-ph]].
  * (36) D. A. Camargo, M. D. Campos, T. B. de Melo and F. S. Queiroz, “A Two Higgs Doublet Model for Dark Matter and Neutrino Masses,” arXiv:1901.05476 [hep-ph].
  * (37) D. Cogollo, R. D. Matheus, T. B. de Melo and F. S. Queiroz, “Type I + II Seesaw in a Two Higgs Doublet Model,” arXiv:1904.07883 [hep-ph].
  * (38) W. C. Huang, H. Ishida, C. T. Lu, Y. L. S. Tsai and T. C. Yuan, “Signals of New Gauge Bosons in Gauged Two Higgs Doublet Model,” Eur. Phys. J. C 78, no. 8, 613 (2018) doi:10.1140/epjc/s10052-018-6067-7 [arXiv:1708.02355 [hep-ph]].
  * (39) W. C. Huang, Y. L. S. Tsai and T. C. Yuan, “Gauged Two Higgs Doublet Model confronts the LHC 750 GeV diphoton anomaly,” Nucl. Phys. B 909, 122 (2016) doi:10.1016/j.nuclphysb.2016.05.002 [arXiv:1512.07268 [hep-ph]].
  * (40) CEPC-SPPC Study Group, “CEPC-SPPC Preliminary Conceptual Design Report. 1. Physics and Detector,” IHEP-CEPC-DR-2015-01, IHEP-TH-2015-01, IHEP-EP-2015-01.
  * (41) P. Q. Hung, “A Model of electroweak-scale right-handed neutrino mass,” Phys. Lett. B 649, 275 (2007) doi:10.1016/j.physletb.2007.03.067 [hep-ph/0612004].
  * (42) H. Ruegg and M. Ruiz-Altaba, “The Stüeckelberg field,” Int. J. Mod. Phys. A 19, 3265 (2004) doi:10.1142/S0217751X04019755 [hep-th/0304245].
  * (43) M. Tanabashi et al. [Particle Data Group], “Review of Particle Physics,” Phys. Rev. D 98, no. 3, 030001 (2018). doi:10.1103/PhysRevD.98.030001
  * (44) B. Kors and P. Nath, “Aspects of the Stüeckelberg extension,” JHEP 0507, 069 (2005) doi:10.1088/1126-6708/2005/07/069 [hep-ph/0503208].
  * (45) B. Kors and P. Nath, “How Stüeckelberg extends the standard model and the MSSM,” doi:10.1142/9789812701756_0056 hep-ph/0411406.
  * (46) B. Kors and P. Nath, “A Supersymmetric Stüeckelberg U(1) extension of the MSSM,” JHEP 0412, 005 (2004) doi:10.1088/1126-6708/2004/12/005 [hep-ph/0406167].
  * (47) B. Kors and P. Nath, “A Stüeckelberg extension of the standard model,” Phys. Lett. B 586, 366 (2004) doi:10.1016/j.physletb.2004.02.051 [hep-ph/0402047].
  * (48) D. Feldman, Z. Liu and P. Nath, “The Stüeckelberg extension and milli weak and milli charged dark matter,” AIP Conf. Proc. 939, 50 (2007) doi:10.1063/1.2803786 [arXiv:0705.2924 [hep-ph]].
  * (49) D. Feldman, Z. Liu and P. Nath, “The Stüeckelberg Z-prime Extension with Kinetic Mixing and Milli-Charged Dark Matter From the Hidden Sector,” Phys. Rev. D 75, 115001 (2007) doi:10.1103/PhysRevD.75.115001 [hep-ph/0702123 [HEP-PH]].
  * (50) D. Feldman, Z. Liu and P. Nath, “The Stüeckelberg $Z$ Prime at the LHC: Discovery Potential, Signature Spaces and Model Discrimination,” JHEP 0611, 007 (2006) doi:10.1088/1126-6708/2006/11/007 [hep-ph/0606294].
  * (51) J. Erler and P. Langacker, “Electroweak model and constraints on new physics,” hep-ph/0407097.
  * (52) S. Schael et al. [ALEPH and DELPHI and L3 and OPAL and SLD Collaborations and LEP Electroweak Working Group and SLD Electroweak Group and SLD Heavy Flavour Group], “Precision electroweak measurements on the $Z$ resonance,” Phys. Rept. 427, 257 (2006) doi:10.1016/j.physrep.2005.12.006 [hep-ex/0509008].
  * (53) U. Baur, O. Brein, W. Hollik, C. Schappacher and D. Wackeroth, “Electroweak radiative corrections to neutral current Drell-Yan processes at hadron colliders,” Phys. Rev. D 65, 033007 (2002) doi:10.1103/PhysRevD.65.033007 [hep-ph/0108274].
  * (54) S. Schael et al. [ALEPH and DELPHI and L3 and OPAL and LEP Electroweak Collaborations], “Electroweak Measurements in Electron-Positron Collisions at W-Boson-Pair Energies at LEP,” Phys. Rept. 532, 119 (2013) doi:10.1016/j.physrep.2013.07.004 [arXiv:1302.3415 [hep-ex]].
  * (55) M. Aaboud et al. [ATLAS Collaboration], “Search for new high-mass phenomena in the dilepton final state using 36 fb-1 of proton-proton collision data at $\sqrt{s}=13$ TeV with the ATLAS detector,” JHEP 1710, 182 (2017) doi:10.1007/JHEP10(2017)182 [arXiv:1707.02424 [hep-ex]].
  * (56) M. Aaboud et al. [ATLAS Collaboration], “Measurement of the cross-section for electroweak production of dijets in association with a $Z$ boson in pp collisions at $\sqrt{s}$ = 13 TeV with the ATLAS detector,” Phys. Lett. B 775, 206 (2017) doi:10.1016/j.physletb.2017.10.040 [arXiv:1709.10264 [hep-ex]].
  * (57) G. Aad et al. [ATLAS Collaboration], “Search for high-mass dilepton resonances using 139 fb-1 of $pp$ collision data collected at $\sqrt{s}=$13 TeV with the ATLAS detector,” arXiv:1903.06248 [hep-ex].
  * (58) J. Alwall et al., “The automated computation of tree-level and next-to-leading order differential cross sections, and their matching to parton shower simulations,” JHEP 1407, 079 (2014) doi:10.1007/JHEP07(2014)079 [arXiv:1405.0301 [hep-ph]].
  * (59) F. Feroz, M. P. Hobson and M. Bridges, “MultiNest: an efficient and robust Bayesian inference tool for cosmology and particle physics,” Mon. Not. Roy. Astron. Soc. 398, 1601 (2009) doi:10.1111/j.1365-2966.2009.14548.x [arXiv:0809.3437 [astro-ph]].
  * (60) K. Fujii et al., “Physics Case for the 250 GeV Stage of the International Linear Collider,” arXiv:1710.07621 [hep-ex].
  * (61) D. d’Enterria, “Physics case of FCC-ee,” Frascati Phys. Ser. 61, 17 (2016) [arXiv:1601.06640 [hep-ex]].
  * (62) S. Martellotti, “The NA62 Experiment at CERN,” arXiv:1510.00172 [physics.ins-det].
  * (63) T. Aushev et al., “Physics at Super B Factory,” arXiv:1002.5012 [hep-ex].
  * (64) S. Alekhin et al., “A facility to Search for Hidden Particles at the CERN SPS: the SHiP physics case,” Rept. Prog. Phys. 79, no. 12, 124201 (2016) doi:10.1088/0034-4885/79/12/124201 [arXiv:1504.04855 [hep-ph]].
  * (65) A. Pilaftsis and C. E. M. Wagner, “Higgs bosons in the minimal supersymmetric standard model with explicit CP violation,” Nucl. Phys. B 553, 3 (1999) doi:10.1016/S0550-3213(99)00261-8 [hep-ph/9902371].
  * (66) V. D. Barger and K. Whisnant, “Heavy $Z$ Boson Decays to Two Bosons in $E(6)$ Superstring Models,” Phys. Rev. D 36 (1987) 3429. doi:10.1103/PhysRevD.36.3429