Consistency of Gauged Two Higgs Doublet Model: Gauge Sector

TL;DR

This paper analyzes constraints on the gauge sector parameters of the gauged two Higgs doublet model using data from LEP, LHC, and projected CEPC sensitivities, clarifying the role of the Stueckelberg mass.

Contribution

It provides the first comprehensive analysis of gauge sector constraints in the gauged two Higgs doublet model using current and future collider data.

Findings

Constraints on new gauge parameters from LEP and LHC data

Projected sensitivities of CEPC on model parameters

Clarification on the zero Stueckelberg mass for hypercharge

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\"{u}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.

Figures

Tables

Table 1. Table 1: Particle content and their quantum number assignments in G2HDM.

Fields	Spin	$SU{\left(3\right)}_C$	$SU{\left(2\right)}_L$	$SU{\left(2\right)}_H$	$U{\left(1\right)}_Y$	$U{\left(1\right)}_X$
$H={\left(H_1{H}_2\right)}^T$	0	1	2	2	$\frac{1}{2}$	$1$
	0	1	1	3	0	0
${\Phi }_H={\left({\Phi }_1{\Phi }_2\right)}^T$	0	1	1	2	0	$1$
$Q_L={\left(u_L{d}_L\right)}^T$	$\frac{1}{2}$	3	2	1	$\frac{1}{6}$	0
$U_R={\left(u_R{u}_R^H\right)}^T$	$\frac{1}{2}$	3	1	2	$\frac{2}{3}$	$1$
$D_R={\left(d_R^H{d}_R\right)}^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)}^T$	$\frac{1}{2}$	1	2	1	$-\frac{1}{2}$	0
$N_R={\left({\nu }_R{\nu }_R^H\right)}^T$	$\frac{1}{2}$	1	1	2	0	$1$
$E_R={\left(e_R^H{e}_R\right)}^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_{\mu }^a\left(a=1,\cdots ,8\right)$	1	8	1	1	0	0
$W_{\mu }^i\left(i=1,2,3\right)$	1	1	3	1	0	0
$W_{\mu }^{\prime i}\left(i=1,2,3\right)$	1	1	1	3	0	0
$B_{\mu }$	1	1	1	1	0	0
$X_{\mu }$	1	1	1	1	0	0

Table 2. Table 2: The electroweak observables at the Z 𝑍 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.

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 }_{\mathrm{had}}$ [GeV]	1.7444 $\pm$ 0.0020	—	1.7411 $\pm$ 0.0008
${\Gamma }_{\mathrm{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 }_{\mathrm{had}}\left[\mathrm{nb}\right]$	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_{\mathrm{FB}}^{(0,e)}$	0.0145 $\pm$ 0.0025	—	0.01618 $\pm$ 0.00006
$A_{\mathrm{FB}}^{(0,\mu )}$	0.0169 $\pm$ 0.0013	—	0.01618 $\pm$ 0.00006
$A_{\mathrm{FB}}^{(0,\tau )}$	0.0188 $\pm$ 0.0017	—	0.01618 $\pm$ 0.00006
$A_{\mathrm{FB}}^{(0,b)}$	0.0992 $\pm$ 0.0016	$0.15\%$	0.1030 $\pm$ 0.0002
$A_{\mathrm{FB}}^{(0,c)}$	0.0707 $\pm$ 0.0035	—	0.0735 $\pm$ 0.0001
$A_{\mathrm{FB}}^{(0,s)}$	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

Equations

$${\mathcal{M}}_{\text{gauge}}^2=\left(\begin{array}{cccc}\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{g{g}_H{v}^2}{4}& -\frac{g{g}_X{v}^2}{2}\\ \frac{g^{\prime }{g}_H{v}^2}{4}& -\frac{g{g}_H{v}^2}{4}& \frac{g_H^2\left(v^2+v_{\Phi }^2\right)}{4}& \frac{g_H{g}_X\left(v^2-v_{\Phi }^2\right)}{2}\\ \frac{g^{\prime }{g}_X{v}^2}{2}+M_X{M}_Y& -\frac{g{g}_X{v}^2}{2}& \frac{g_H{g}_X\left(v^2-v_{\Phi }^2\right)}{2}& g_X^2\left(v^2+v_{\Phi }^2\right)+M_X^2\end{array}\right).$$

$$({\mathcal{O}}^{4\times 4}{)}^{\mathrm{T}}\cdot {\mathcal{M}}_{\text{gauge}}^2\cdot {\mathcal{O}}^{4\times 4}=\text{diag}(0,M_{Z_1}^2,M_{Z_2}^2,M_{Z_3}^2),$$

$${\mathcal{O}}_{M_Y=0}^{4\times 4}=\left(\begin{array}{cccc}{c}_W& -s_W& 0& 0\\ s_W& c_W& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right)\cdot \left(\begin{array}{cccc}1& 0& 0& 0\\ 0& & & \\ 0& & \mathcal{O}& \\ 0& & & \end{array}\right),$$

$$e^{\mathrm{SM}}\equiv g\sin {\theta }_W=g^{\prime }\cos {\theta }_W=\frac{g{g}^{\prime }}{\sqrt{g^2+g^{\prime 2}}}.$$

$${\mathcal{O}}_{\text{SM}}^{4\times 4\mathrm{T}}\cdot {\mathcal{M}}_{\text{gauge}}^2(M_Y=0)\cdot {\mathcal{O}}_{\text{SM}}^{4\times 4}$$

$$\mathcal{O}=\left(\begin{array}{ccc}{c}_{\psi }{c}_{\varphi }-s_{\theta }{s}_{\varphi }{s}_{\psi }& -c_{\theta }{s}_{\varphi }& s_{\psi }{c}_{\varphi }+s_{\theta }{s}_{\varphi }{c}_{\psi }\\ c_{\psi }{s}_{\varphi }+s_{\theta }{c}_{\varphi }{s}_{\psi }& c_{\theta }{c}_{\varphi }& s_{\psi }{s}_{\varphi }-s_{\theta }{c}_{\varphi }{c}_{\psi }\\ -c_{\theta }{s}_{\psi }& s_{\theta }& c_{\theta }{c}_{\psi }\end{array}\right),$$

$$\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{)}}\,,$$

$$\tan (\theta )=\frac{-g_X(M_{Z_2}^2(v^2-v_{\Phi }^2)+M_{Z^{\text{SM}}}^2{v}_{\Phi }^2)}{v{M}_{Z^{\text{SM}}}(M_X^2-M_{Z_2}^2+2g_X^2{v}_{\Phi }^2)}\sin \varphi ,$$

$$\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 \varphi }-\sin \theta \cot \varphi .$$

$$(Z_1,Z_2,Z_3{)}^{\mathrm{T}}={\mathcal{O}}^{\mathrm{T}}\cdot (Z^{\mathrm{SM}},W^{\prime 3},X{)}^{\mathrm{T}}.$$

$${\mathcal{L}}_N=g_M\sum _f\sum _{i=1}^3\bar{f}{\gamma }_{\mu }\left[\left(v_f^{(i)}-{\gamma }_5{a}_f^{(i)}\right)Z_i^{\mu }\right]f,$$

$$v_f^{(i)}=$$

$$+\frac{1}{\sqrt{g^2+g^{\prime 2}}}\left(X_R{g}_X{\mathcal{O}}_{4,i+1}^{4\times 4}+T_{f_R}^{3H}{g}_H{\mathcal{O}}_{3,i+1}^{4\times 4}\right),$$

$$a_f^{(i)}=$$

$$-\frac{1}{\sqrt{g^2+g^{\prime 2}}}\left(X_R{g}_X{\mathcal{O}}_{4,i+1}^{4\times 4}+T_{f_R}^{3H}{g}_H{\mathcal{O}}_{3,i+1}^{4\times 4}\right).$$

$${\mathcal{L}}_{\gamma }=-e^{\text{SM}}\sum _f\bar{f}{\gamma }_{\mu }\left(Q_f^{G2HDM}-a_f^{\gamma }{\gamma }_5\right)A^{\mu }f,$$

$$Q_f^{G2HDM}=$$

$$+\frac{1}{2e^{\text{SM}}}\left(g_X{\mathcal{O}}_{4,1}^{4\times 4}{X}_R+g_H{\mathcal{O}}_{3,1}^{4\times 4}{T}_{f_R}^{3H}\right),$$

$$a_f^{\gamma }=$$

$$-\frac{1}{2e^{\text{SM}}}\left(g_X{X}_R{\mathcal{O}}_{4,1}^{4\times 4}+g_H{T}_{f_R}^{3H}{\mathcal{O}}_{3,1}^{4\times 4}\right).$$

$$v_{f(M_Y=0)}^{(i)}=$$

$$a_{f(M_Y=0)}^{(i)}=$$

$$v_f^i=$$

$$a_f^i=$$

$$\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{]}\,,$$

$${\mathcal{R}}_f=\left(1+{\delta }_f^{QED}\right)\left(1+\frac{N_f^c-1}{2}{\delta }_f^{QCD}\right),$$

$${\delta }_f^{QED}$$

$${\delta }_f^{QCD}$$

$$R_l=\frac{{\Gamma }_{\mathrm{had}}}{{\Gamma }_{l^{+}{l}^{-}}},R_q=\frac{{\Gamma }_{q\bar{q}}}{{\Gamma }_{\mathrm{had}}},$$

$${\sigma }_{\mathrm{had}}=\frac{12\pi {\Gamma }_{e^{+}{e}^{-}}{\Gamma }_{\mathrm{had}}}{M_Z^2{\Gamma }_Z^2},$$

$$A_f=\frac{2v_f{a}_f}{v_f^2+a_f^2},$$

$$A_{\mathrm{FB}}=\frac{3}{4}{A}_f\frac{A_e+P_e}{1+P_e{A}_e},$$

$$A_{\mathrm{FB}}^{(0,f)}=\frac{3}{4}{A}_e{A}_f.$$

$${\chi }_{Z-\mathrm{pole}}^2$$

$${\chi }_i^2=\frac{{\left(O_i^{\exp }-O_i^{\mathrm{th}}\right)}^2}{{\left(\Delta O_i^{\exp }\right)}^2},$$

Full text

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}$.

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

[TABLE]

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}$

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

[TABLE]

[TABLE]

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

[TABLE]

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

[TABLE]

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$)

[TABLE]

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

[TABLE]

where

[TABLE]

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

[TABLE]

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)

[TABLE]

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

[TABLE]

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

[TABLE]

with

[TABLE]

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

[TABLE]

the hadronic cross-section

[TABLE]

the parity violation quantity

[TABLE]

and the forward-backward asymmetry quantity

[TABLE]

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

[TABLE]

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

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

[TABLE]

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

[TABLE]

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

[TABLE]

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:

[TABLE]

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

[TABLE]

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.

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

[TABLE]

where

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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,

[TABLE]

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.,

[TABLE]

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

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.

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.

IV.3 Light $M_{X}$ Scenario

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.

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}$.

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).

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).

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

[TABLE]

where

[TABLE]

and

[TABLE]

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.

[TABLE]

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.

[TABLE]

[TABLE]

[TABLE]

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

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

• •

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

[TABLE]

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

[TABLE]

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

[TABLE]

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}$.