Consistency of Gauged Two Higgs Doublet Model: Gauge Sector

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.


I. INTRODUCTION
The discovery of the 125 GeV scalar boson identified as the Higgs boson in the Standard Model (SM) [1][2][3][4] suggested that the simple Higgs mechanism [5][6][7] for electroweak symmetry breaking proposed by Weinberg [3] and Salam [4] 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) [8] and dark energy [9]. The standard ΛCDM model of cosmology [10] 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 [11] was motivated partly by the inert Higgs doublet model (IHDM) [12][13][14][15] of dark matter. IHDM is a variant of the general 2HDM [16] with an imposed discrete 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 [17][18][19][20][21][22][23][24][25][26]. In G2HDM, the discrete 2 symmetry in IHDM was not enforced. Instead the two Higgs doublets 1 and 2 are grouped into a twodimensional irreducible representation = ( 1 , 2 ) T of a new gauge group (2) .
A priori there is no need to impose the discrete 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 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 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 [27,28]. Detailed analysis of the complex scalar dark matter physics in G2HDM will be presented in a forthcoming paper [29].
The construction of G2HDM in [11] 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 [30,31,38,39]. In particular, we have studied recently in details the theoretical and phenomenological constraints on the scalar sector [31].
We note that the 2HDM augmented with an extra local abelian (1) has been discussed in the literature [32][33][34][35][36][37] 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 [31] done for the scalar sector. It is also interesting to address the sensitivities of these new parameters at the future colliders. In this section, we will start with a brief review for the set-up of G2HDM [11] by specifying its particle content (Sec. II A) and then write down the mass spectrum of the neutral gauge bosons (Sec. II B) and their interactions with the SM fermions (Sec. II C) in the model. Along the way, we will discuss some peculiar effects for nonzero Stüeckelberg mass associated with the hypercharge (1) .

A. Particle Content
The particle content of G2HDM is listed in in G2HDM, we refer our readers to [11] 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 (1) charges. Thus the complete gauge groups in G2HDM consist of Here , ′ , and denote the gauge couplings of  (1) and (1) respectively. We note that Δ the VEV of the triplet Δ does not enter into the neutral gauge boson mass matrix. This is unlike the case of scalar boson mass matrix analyzed in [31] which involves all three VEVs, , and Δ . The matrix ℳ 2 gauge in Eq. (1) is real and symmetric and thus can be diagonalized by a 4×4 orthogonal rotation matrix that we will denote as 4×4  [43], which implies must be very tiny too. If individual Stüeckelberg mechanism is introduced for each of the two (1)s factors in G2HDM, the photon will in general obtain nonzero mass and many results obtained in [42] apply as well.
In [11], we followed [44][45][46][47] in which only one Stüeckelberg field was introduced for the two factors of (1)s to implement the Stüeckelberg mechanism. The matrix ℳ 2 gauge 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 and .
In the next subsection, we will show that with a nonzero 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 [42].
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 through the mixing matrix elements and hence necessarily small. Thus, we have strong theoretical motivation to set = 0 in what follows to avoid these unpleasant features. For an analysis with both and nonzero in a Stüeckelberg (1) extension of the SM that maintains the standard QED interaction for the SM fermions, see [48][49][50]. 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 (2) whereas there is only one extra abelian group (1) in [48][49][50].
Setting = 0 in G2HDM will simplify ℳ 2 gauge and allows us to write the rotation matrix in the following product form where and represent cos and sin respectively, with being the Weinberg angle defined by It is obvious that the matrix 4×4 =0 in Eq. (3) is just the product of the SM gauge rotation matrix made into a 4×4 matrix, called 4×4 SM , times a general 3×3 orthogonal rotation matrix which was also converted to a 4 × 4 matrix. After applying the where SM = √︀ 2 + ′2 /2 is the mass of the boson in the SM. We can consider the vanishing (1,1) element to be the mass of the photon eigenstate . Furthermore, according to Eqs. (2) and (3), the remaining 3×3 matrix formed by the non-vanishing elements above is diagonalized by the orthogonal matrix . In particular, one can parametrize in terms of the following Tait-Bryan representation where and stand for sine and cosine with the rotation angle = , , respectively. As shown in Appendix A, these rotation angles can be represented as It is easy to see that taking the limits of and go to 0, the non-vanishing 3×3 block matrix in Eq. (5) becomes Diag( 2 SM , 0, 2 ). Thus the rotation matrix must be identity. This can be realized by setting , and to be zeros which can be derived from Eqs. (7), (8) and (9).
We note that if one sets to zero, the mass matrix in the right-handed side of Eq. (5) is symmetric under the interchange of /2 ↔ .
After the rotation matrix is found, the mass eigenstates where runs from 1 to 3 are given by The composition SM , ′3 and of the mass eigenstate is given by 2 1 , 2 2 , and 2 3 , respectively. In general, the -pole can be any one of the 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 -pole.

C. Neutral Gauge Current Interactions
The part of the Lagrangian that contains the interaction of the with visible matter in G2HDM is where = √︀ 2 + ′2 /2. The ( ) and ( ) factors are given by ( Here 3 is the For the photon-fermion couplings in G2HDM, we obtain where Thus, with both nonzero and , 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 4×4 1,1 / plus two correction terms, and there is also a non-vanishing axial vector coupling. Next, we can take the limit = 0 and write the corresponding expressions. By replacing the elements of 4×4 by 4×4 1,1 = 4×4 2,2 = and − 4×4 1,2 = 4×4 2,1 = as in Eq. (3), one can find the following new expressions for the vector and axial vector Similarly, one can do the same substitutions on Eqs. (15) and (16)  vanish. This is the main physical reason why we set = 0 so as to reproduce the standard photon-fermion couplings. For , it can be arbitrary and is naturally to consider the light and heavy scenarios where it is smaller and greater than the -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.

A. Constraints from Precision Electroweak Data at LEP-I
The interaction of boson with SM fermions is described by the Lagrangian in Eq. (11). For the case of = 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 boson and fermions couplings [51,52], which now are given by 2 (suppressing = 0 in the subscripts) where in this work is either equal to 1 or 2 depending which mass eigenstate is closest to -pole. The parameters and are loop corrections quantities. The decay of the boson into fermions and anti-fermions in the on-shell renormalization scheme is given by [51,53] where is the color factor (1 for leptons and 3 for quarks), Γ = 3 /6 √ 2 , Here is the electric charge of the fermion in unit of SM , and and are the fine-structure and strong coupling constants, respectively, evaluated at the scale.
It is understood that the couplings and in Eq. (21) should be replaced by and in Eqs. (19) and (20) respectively with = 1 or 2 depending which is closest to the -pole .
We also investigate some -pole ( √ ≈ ) observables, including the ratio of partial decay width of boson the hadronic cross-section the parity violation quantity and the forward-backward asymmetry quantity where is the initial − polarization. Recall that at LEP-I = 0, in this case is presented in Table II.
From the data in Table II, we build the Chi-squared for the electroweak observables at -pole as follows Note that we have considered the correlations between the total decay width of boson and its partial decay widths to hadrons, invisibles and dilepton. For each 2 on the right-handed side of Eq. (30), it is given by the standard expression, namely where exp/th represents the experimental/theoretical value of any one of the 23 electroweak observables listed in Table II and Δ exp is the corresponding experimental uncertainty.

B. Contact Interactions at LEP-II
We also include constraints from data above the -pole by considering the LEP-II measurements related to contact interactions taking the following form of effective where , represent the chirality projection operators with , being or 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  [43], CEPC preliminary conceptual design report [40], and the SM prediction [43], respectively.
for the indices of Λ ± : = { , , , , , }, which are also called models. The limits on Λ ± set by LEP-II are given in Table 3.15 of Ref. [54]. The strongest constraint is given by Λ + > 24.6 TeV. By using these Λ ± 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 -pole, we calculate the additional -like mediator contribution 3 to the − + → →¯scattering cross section. In the case = we have the contribution of both and channels while for ̸ = only the channel contributes. Note that here we do not need the SM contributions such as the photon and exchange not considered in Eq. (32). In the massless approximation for all the external fermions, the amplitudes for the and channels and for the interference term between them are given by: where is the center of mass energy squared, = (cos − 1)/2 and is the angle between incoming and outgoing particles. This angle 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 Λ ± given by the experimental result. The couplings and have = 1 or 2 depending on whether we are analyzing light or heavy scenario. For = 3, we assume 3 is much heavier than so that its contributions are negligible. To be able to construct a 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 2 value should be 2.71 4 . In this case, we calculate the 2 value using where eff is the cross section obtained using the effective Lagrangian of Eq. (32) with the experimental results for Λ ± given in Ref. [54] for  distinguish signal from background. It is commonly believed that the Drell-Yan con-straint¯→ → + − 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 [55]. The data of electron-positron pair ( ) and muon-pair ( ) final states are given by Tables 3 and 4 respectively in Ref. [55]. In the signal region located around -boson mass (the invariant mass 80 < / GeV < 120), we found that the systematic uncertainties of Drell-Yan background is larger than the data statistic uncertainties in both or final state. We have also checked that the EWPT constraints in Table II are much stronger than LHC Drell-Yan constraint. On the other hand, -boson can be singly produced either by radiation from the incoming partons (Fig. 1a) or -channel exchange of a gauge boson (Fig. 1b).
To constrain the G2HDM modified + − couplings, the later process is more use-ful than the former because QCD processes usually suffer from larger systematical uncertainties than the electroweak ones. Recently, ATLAS [56] reported a fiducial electroweak cross section of EW = 119±16±20±2 fb and EW = 34.2±5.8±5.5±0.7 fb for dijet invariant masses greater than 250 GeV and 1 TeV, respectively. The SM simulated cross sections EW (SM) are also given in Table 5  Comparing with the SM, except for the + − 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 where and = , . However, similar to direct Drell-Yan boson search, we found that the value of ℛ 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 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 decay in our numerical works so as to save some computer resources.

LHC ′ 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 [38]. Unlike the study in Ref. [38] which In Fig. 3 of Ref. [57], one can see the upper limits of cross section times branching ratio BR( ′ → + − ) are based on the ratio of the total width Γ ′ of ′ divided by its ′ . Depending on this ratio, the limits can be altered by a factor of ∼ 5. As shown in Appendix B, the Γ / 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 Γ ′ / ′ = 0.06. 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 = /10 for maximum invisible decay but we found that the Γ( → * ) can differ within an accepted range of ∼ 6% comparing with the massless case.
Using  (17) and (18). Hence, by using the same reasoning as before we include the latest ATLAS ′ limit in our scan by using the following chi-squared function where the branching ratio BR G2HDM ( → + − ) can be found in Appendix B and 95% ATLAS ( → ′ )×BR 95% ATLAS ( ′ → + − ) is 95% C.L. taken from the curve associated with Γ ′ / ′ = 0.06 in Fig. 3 of Ref. [57].

A. 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 where is the smallest area bound with a fraction of the total probability and the normalization in the denominator is the total probability with → ∞.
The total 2 Total ( , , Φ , ) we will use in our analysis is the sum of Eqs.
where we have suppressed the arguments of all the 2 functions. We adopt the statistical sensitivity as Since our likelihood is modeled as a pure Gaussian distribution, i.e. ℒ ∝ exp(− 2 /2), one can connect the 2 to the confidence level: the 68% (95%) C.L. in a two dimensional parameter space corresponding to Δ 2 = −2 ln(ℒ/ℒ max ) = 2.30 (5.99). Here ℒ max is the maximum value of the likelihood in the region .
There are two interesting scenarios: (i) heavy and (ii) light . The heavy scenario will result in two new heavy neutral gauge bosons 2 ≡ ′ and 3 ≡ ′′ , and the measured boson located at -pole will be the lightest one, 1 ≡ . However, the light scenario will result in a new boson 1 lighter than the -pole which is usually called dark ( ) or dark photon ( ). In this case, 2 corresponds to the -pole 2 ≡ and 3 ≡ ′ . Hence, we choose our scan ranges for two scenarios, .
For the other three parameters, we use the same ranges for the two scenarios of 5 , We perform random scans by using MultiNest v2.17 [59] 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 ∼ 10 5 samples for each scenario.

B. Heavy Scenario
In the heavy scenario, the mass of 1 boson is located at around -pole (∼ 91 GeV) so that 1 is identified as the SM -boson. Note that 1 ( ) boson physics is strongly affected by the different composition of 2 ( ′ ) but not the heaviest boson 3 ( ′′ ) because 3 is heavier than 2 in our parameter choices and therefore has less impact.
In Fig. 2, we present the scatter points of the composition of 2 = 12 SM + 22 ′3 + 32 for the 1 region based on the likelihoods described in Sec. IV A.
The color code hereafter represents the three different composition of 2 . Recalling Eq. (10), we define ′3 -like 2 with condition 2 22 > 0.8 (red crosses ×), mixed state The 1 allowed scatter points projected on the ( 2 , 2 22 ) and ( 2 , 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 2 (blue triangles) is less evenly distributed because it needs some trade-off between the two new gauge couplings and . In Fig. 2b, we projected the same parameter space on the plane ( 2 , 2 12 ). Note that the mixing 2 12 presents how 2 is consisted of SM . Therefore, very small 2 12 implies 2 32 ≈ (1 − 2 22 ) from the orthogonality of . Furthermore, the upper limit of Φ sets an lower limit of the 2 12 for the red cross region. If Φ goes to infinity, 2 becomes super heavy and decouple. The composition of SM in 2 should then be negligible, thus 2 12 vanishes in this limit. We note that the excluded concave up region of 2 between 250 GeV and 6 TeV on the upper limit of 2 12 is due to the constraint from ATLAS ′ search.
In Fig. 3, we show the 1 (dashed) and 2 (solid) likelihood contours with scatter points inside the 1 region on the (a) ( 2 , ) and (b) ( 2 , Φ ) planes. In Fig. 3a, we can see that the ′3 -like red crosses form a band with a tendency proportional to . This is because for a ′3 -like 2 , 2 2 ≈ 2 ( 2 + 2 Φ )/4 ≈ 2 2 Φ /4 which can be extracted from the (3,3) element of the mass matrix in Eq. see the green circles running from below the two red cross and blue triangle bands up to the upper limit of . In other words, we can see how the 2 passes from being dominated by the (3,3) element of Eq. (5) (red crosses), which is -dependent, to being dominated by the -independent (4,4) element (green circles).
One particular feature of Fig. 3b is that the low 2 and low Φ region (lower left corner) is covered only with -like points while both ′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 ′3 and mixed points in Fig. 3a for low 2 < 200 GeV and 10 −2 . These curves in the upper bound (Fig. 3a) and in the lower bound (Fig. 3b)  is not displayed for the -like points since they do not depend strongly on .
The ATLAS ′ constraint almost rules out the region 250 GeV < 2 < 6 TeV for ′3 -like and mixed 2 , except the region with Φ > 100 TeV. However, the -like 2 at the same region has not been affected much by the ATLAS ′ constraint.
Similarly, in Fig. 4, we show the 1 (dashed) and 2 (solid) likelihood contours with scatter points inside the 1 region on the (a) ( 3 , ) and (b) ( 3 , ) planes. From Fig. 4a, one can easily see that the -like 2 boson (green circles) forms a band whose tendency is proportional to the . This can be understood by the fact that the composition of the 3 in this case is mainly from ′3 , which has a mass proportional to 0.5 √︀ 2 + 2 Φ ≈ 0.5 Φ again coming from the (3,3) element of the mass matrix in Eq. (5). On the other hand, in the case of the ′3 -like 2 boson (red crosses), the mass of the 3 almost does not depend on . Indeed, the composition of 3 is now mainly from and 2 3 ≈ ( 2 ( 2 + 2 Φ ) + 2 ). This is clearly shown in Fig. 4b, when is small ( < 3 × 10 −3 ), the mass of 3 in the red cross region is dominated by and less than our set-up limit of 10 4 GeV. However, when is getting bigger, the mass of the 3 can be dominated by the Φ term for sufficiently large value of Φ . We can also see that the EWPT data sets upper bounds on and . The excluded concave up region of 250 GeV < 2 < 6 TeV in Fig. 4a for the ′3 -like and mixed composition of 2 is again due to the ATLAS ′ search which does not apply for the -like case. As a result, the ATLAS ′ search cannot constrain on for ′3 -like points as clearly shown in Fig. 4b.

C. Light Scenario
In the light scenario, we require that the mass of 2 boson is always at around -pole (∼ 91 GeV). In this scenario, the lightest 1 with mass less than the -boson mass can be the dark photon or dark , while the conventional ′ is the heaviest boson 3 . We note that the composition of 3 is given by 3 = 13 SM + 23 ′3 + 33 . The 1 allowed scatter points projected on the ( 3 , 2 23 ) and ( 3 , 2 13 ) planes are depicted in Figs. 5a and 5b, respectively. The color code for the composition of 3 is the same as in Fig. 2 for 2 .
An obvious feature of Fig. 5a is that the mixed state of 3 (blue triangles) has a mass upper limit. Intuitively, it requires some trade-off between the gauge couplings GeV. This effect will be discussed with more detail later in Fig. 6. In Fig. 5b, we can see that the SM composition of 3 is again small. However, unlike the heavy scenario, the -like 3 boson has a similar distribution as ′3 -like 3 boson. Additionally, the mixed 3 state at the mass region between 210 GeV and 700 GeV cannot be excluded by the ATLAS ′ constraint which is also different from the heavy scenario.
In analogous to Fig. 3, we show in Fig. 6   ) planes. Again, the red cross represents the points of ′3 -like 3 boson, the blue triangle represents the points of mixed state ( SM , ′3 and ) 3 boson, and the green circle represents -like 3 boson. We note that in this scenario, 1 is considered as a dark photon 6 and has mass range from 1 MeV to -pole. One can easily see that the 3 composition is clearly separated on the planes of ( 1 , ) and

D. 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 -boson precision experiments can further probe our model. In the near future, there are three colliders that can improve -boson measurements: CEPC [40], ILC [60], and FCC-ee [61]. 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 6 The lightest boson 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. the projected CEPC sensitivity.
In the third column of Table II, we quote the expected CEPC sensitivity [40].
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 Δ 2 in terms of Φ in Fig. 8  In Fig. 9, we compare the present limit and future CEPC sensitivity of the twodimensional contours on the ( , ) plane. The figure in the left (right) column corresponds to the heavy (light) scenario. Because the upper scan limit of is set to be less than SM 2 , the experimental constraints on are not present. In contrast, has an upper limit from the constraints due to having a lower limit from the maximum scanned Φ value. The upper limit of can be further improved by future CEPC sensitivity along the edge of the contour. However, the light 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 ′3like 2 in heavy scenario (left panel) or ′3 -like 3 in light scenario (right panel), has a lower limit at ∼ 2 × 10 −3 due to our choices of the parameter scan ranges. Indeed, in both cases we have 2,3 ≈ 0.5  it is again too light to be presented in the EWPT data. To constrain light , just like dark photon, the lightest 1 could be detected by those future beam dump experiments such as NA62 [62], Belle II [63], and SHiP [64]. 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 and associated with the hypercharge (1) and the extra (1) respectively, we showed that a nonzero would produce non-standard QED couplings for all the fermions in G2HDM, albeit we can always achieve a massless photon for arbitrary values of and . We therefore set = 0 in our numerical analysis. The remaining new parameters in the gauge sector of G2HDM needed to be constrained are , , Φ and .
We have examined the remaining parameter space with the EWPT LEP data at the -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 Φ , 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 decay does not impose noticeable impacts yet.
We classify our parameter space based on three different composition ( -like, ′3 -like, and mixed) of the heavy neutral gauge boson, either 2 or 3 , which is the next-heavier boson than the SM one, in order to manifest the physics and constraints discussed in this paper.
In the heavy scenario ( > 100 GeV), the SM-like is the lightest 1 boson and EWPT constraints exclude the small Φ region up to 24 TeV at 2 significance.
However, the EWPT constraints are not so sensitive to the light scenario ( < 80 GeV) where SM-like is the next-lightest 2 boson. In particular, the Φ is required to be greater than 15 TeV due to the constraints of ′ contact interaction search from LEP-II and high-mass dilepton resonance search from LHC Run II.
Furthermore, in both light and heavy scenarios, 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 in the light scenario. On the other hand, it is not so trivial for the couplings and because we found it is hard to set an upper bound on them individually.
Although the SM boson is fixed at the -pole, the allowed physical masses of the heavier still depend on the and detailed composition. Generally speaking, the 2 allowed mass range in the heavy scenario is same as the range of have no difference between different composition. However, regarding to 3 , mixed 3 is restricted to less than 500 GeV but the masses of -like and ′3 -like 3 are below 70 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 scenario but the sensitivity is not improved much for the light scenario.
In the latter case, when is getting very light, 1 can be much lighter than the -boson and it is more appropriate to identify it as the dark photon or dark .
The very rich phenomenology of light dark photon or dark in G2HDM remains to be explored in the future. 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 ℳ 2 gauge ( = 0) given in Eq. (5). The orthogonal matrix we choose is Eq. (6) because it is rather convenient to find all the s and determine the rotational angles , and 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 , so we apply Cramer's rule for solving the secular equations and get another form for as follows where with the range for covers radians, and the range for and covers 2 radians.
Note that the expressions in Eq. (A4) do not depend on the Δ given in Eq. (A2).
• The new neutral gauge boson can also decay into pair of scalar dark matter candidate in this model. The decay width for this process → * is given by [66] Γ( → * ) = 2

48
( where the coupling = 2 and = 2 2 . We note that is a tripletlike 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 → + − is given by where ± = 2 ± 2 and the coupling + − is given as follows • The decay width for → is given by [66] Γ( → ) = and the coupling is given as follows here is the VEV of the SM Higgs field. Note that we have ignored the mixing of SM like-Higgs with other scalar bosons in the above calculations.
• Finally, if not kinematically prohibited, the new neutral gauge bosons can also decay into ′ and the dark matter . The decay width for this process can be computed as where the coupling ′ = 2 2 Δ with Δ being the VEV of (2) 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 and light scenarios respectively.
In those plots, we set the dark matter mass to be 10% of the new heavy neutral gauge boson ( . . Moreover, we assume the masses of new heavy fermions are degenerate and equal to 3 TeV. We note that Δ can be derived from other parameters according to Δ = 0.5 . From these scatter plots one can see that for the heavy neutral gauge bosons in both scenarios, their ratios Γ / are all below ∼ 1%, until they are heavier than 10 TeV the ratios can then reach ∼ 6%. However for the light 1 in the light scenario, Γ 1 / 1 is well below 10 −4 .