Collider Bounds on 2-Higgs Doublet Models with $U(1)_X$ Gauge Symmetries

2-Higgs Doublet Models (2HDMs) typically need to invoke an ad-hoc discrete symmetry to avoid severe flavor bounds and in addition feature massless neutrinos, thus falling short of naturally complying with existing data. However, when augmented by an Abelian gauge symmetry naturally incorporating neutrino masses via a type-I seesaw mechanism while at the same time escaping flavor changing interactions, such enlarged 2HDMs become very attractive phenomenologically. In such frameworks, the distinctive element is the $Z'$ gauge boson generated by the spontaneous breaking of the Abelian group $U(1)_X$. In this work, we derive updated collider bounds on it. Several theoretical setups are possible, each with different and sometimes suppressed couplings to quarks and leptons. Thus, complementary data from dijet and dilepton resonance searches need to be considered to fully probe these objects. We employ the corresponding datasets as obtained at the Large Hadron Collider (LHC) at the 13 TeV CMs energy for $\mathcal{L}=12,36$ and $300$ fb$^{-1}$ of luminosity. Moreover, we present the potential sensitivity to such $Z'$s of the High Luminosity LHC (HL-LHC) and High Energy LHC (HE-LHC).


I. INTRODUCTION
The Standard Model (SM) offers the best description of strong and Electro-Weak (EW) interactions and has successfully passed all precision tests up to now [1]. In the SM, fermion masses are generated via the presence of one scalar doublet which gives rise to the 125 GeV Higgs boson discovered some years ago at the LHC [2,3]. One of the precision measurements that attest the predictability of the SM concerns the ρ parameter, ρ = M 2 W /(M 2 Z cos θ 2 W ), where θ W is the Weinberg angle. Today, the EW tests lead to ρ = 1.01032 ± 0.00009 [4], with the error bar being driven mostly by the uncertainty of the top quark mass which appears at one loop level in the calculations.
In general extended scalar sectors the ρ parameter is given by where I i and Y i are the isospin and hypercharge of a scalar representation with Vacuum Expectation Value (VEV) v i . Hence models with extended Higgs sectors that feature scalar doublets with hypercharge equal to unity or scalar singlets with zero hypercharge straightforwardly preserve ρ = 1 at tree level. From this perspective, models with the presence of a second Higgs doublet stand out because they also offer a hospitable environment for EW phase transition * dacamargov@gmail.com † luigi.dellerose@fi.infn.it ‡ s.moretti@soton.ac.uk § farinaldo.queiroz@iip.ufrn.br [5][6][7][8][9], collider [10] and flavor physics [11][12][13] (see [14] for a review).
However, such 2HDMs [15] suffer from severe flavor bounds because the presence of a second Higgs doublet induces flavor changing neutral currents (FCNCs) [16,17]. One way around this problem has been to work with Yukawa couplings that obey certain relations so that FCNCs are suppressed [18][19][20]. An orthogonal solution to this problem has been to enforce an ad-hoc Z 2 symmetry under which one scalar doublet is even and the other is odd. That said, it would be elegant if one could realize this discrete symmetry via gauge principles and, in addition, also accommodate neutrino masses since their explanation is still absent in the original formulation of the 2HDM [21][22][23][24][25][26].
Having that in mind, the relation between continuous symmetries and discrete symmetries has been addressed in [18,27]. Gauge symmetries have been considered in several 2HDM scenarios in [28][29][30][31][32], but the aforementioned flavor problem was addressed in this context only in [33] and was later expanded in [34] with additional U (1) X models. In the presence of an additional U (1) X group, neutrino masses are naturally taken into account since the extra gauge symmetry requires three right-handed neutrinos in order to cancel the gauge anomalies [33,34]. Indeed, the corresponding Majorana mass term leads, through a type-I seesaw mechanism, to a simple explanation of the neutrino mass problem. 2HDMs with a U (1) X gauge symmetry are also characterized by a richer phenomenology due to the presence of a massive Z gauge boson that arises after the spontaneous symmetry breaking. Moreover, neutrino masses and dark matter have been addressed recently with gauge symmetries and a singlet scalar extension in [35]. Our goal in this work is to use dijet [36] and dilepton [2,[37][38][39] data from the LHC to constrain the mass of such Z boson and, consequently, the viable parameter space of the model. Instead of focusing on only one U (1) X realization, we will investigate all eight U (1) X models introduced in [34] which are capable of solving the flavor problem in the 2HDM as well as generating neutrino masses from gauge principles.
Our work is structured as follows. In Section II we review the models, in Section III we derive the aforementioned collider limits, in Section IV we comment on the sensitivity at the HL/HE-LHC. We draw our conclusions in Section V.

II. THE MODEL
The 2HDM is characterized by a rich collider and flavor phenomenology but suffers, in general, from FCNC bounds and lacks an explanation for neutrino masses [15]. The 2HDM embedding Abelian gauge groups that solve these problems appeared in [32][33][34]. There are several possible U (1) X symmetries that can be incorporated in the 2HDM that suppress FCNCs (see table I) while being free from triangle anomalies through the addition of three right-handed neutrinos. If the U (1) X charges of the two Higgs doublets are different, the Abelian gauge symmetry naturally prohibits one of the two doublets to participate in the generation of the SM fermion masses. This mechanism replaces the ad-hoc Z 2 discrete symmetry, which is commonly employed in the 2HDM, and provides a Yukawa structure similar to the type I scenario. Furthermore, one can properly choose the U (1) X quantum numbers of all particles such that no anomalies are present. In this respect, three right-handed neutrinos are necessary, which charges are fixed to −(u + 2d), where u and d are the quantum numbers of the right-handed up-and down-quarks, respectively. For instance, if u = d = 1 3 , we end up with the well-known B −L gauge symmetry. However, as seen in table I, many other models are also possible. The SM fermions and the neutrinos acquire mass through the Lagrangian Notice that only the scalar doublet Φ 2 is relevant for the mass generation of SM fermions while Φ 1 is decoupled from the latter. The doublets can be written as, with v 2 = v 2 1 + v 2 2 being the EW VEV. The singlet scalar Φ s = 1/ √ 2 (v S + ρ s + iη s ) is paramount to neutrino masses which are dynamically generated via the last two terms of Eq.
(2) by the singlet VEV v S . It leads to a mass matrix typical of the usual type I seesaw mechanism [40][41][42], which results in m ν = −m T D /M R m D and m N = M R , as long as In order to generate active neutrino masses at the sub-eV scale one can either set the right-handed neutrino masses at the scale of a Grand Unified Theory (GUT) or else adopt suppressed Yukawa couplings [43,44]. The right-handed neutrinos, being charged under U (1) X , may have a relevant impact on the collider bounds of the Z because the latter can decay into them [45][46][47][48]. In this work we will conservatively assume all right-handed neutrinos to have the same mass of 100 GeV. This assumption is conservative in the sense that, if the right-handed neutrinos were sufficiently heavy to prohibit the Z boson to decaying into them, the Branching Ratio (BR) of the Z into charged leptons or light quarks would be larger, thus strengthening our limits.
The scalar fields of the model are described by the following potential Notice the absence in the potential, due to U (1) X invariance, of the m 2 12 Φ † 1 Φ 2 quadratic term. In the standard Z 2 realization of the 2HDMs, one has to introduce an ad-hoc m 2 12 parameter that softly breaks the discrete symmetry. In these scenarios, instead, it is dynamically generated by the vev v S of the singlet scalar. The presence of a scalar doublet charged under U (1) X (see table I) leads to Z − Z mass mixing. This mass mixing is explicitly derived in Appendix A. Moreover, we also account for the presence of a kinetic mixing in the Lagrangian [49][50][51][52][53]: Here, B and X are the neutral vector bosons from the U (1) Y and U (1) X gauge groups, which, after EW Symmetry Breaking (EWSB) and together with the third component of the SU (2) L gauge bosons, give rise to the massless photon as well as massive Z and Z .
In summary, after taking into account the kinetic and mass mixings, we find the neutral current where ξ is the mixing parameter explicitly given in Eq. (29) while Q L X (Q R X ) are the left-handed (right-handed) fermion charges under U (1) X defined according to table I. One can then easily obtain the Z interactions with the SM fermions by substituting their charges for each of the models exhibited in table I. This interaction Lagrangian represents the key information for the collider phenomenology we are going to tackle, because it dictates Z production as well as its most prominent decays to be searched for. In principle, there are other interactions involving the Z gauge boson besides the neutral current of Eq. 9. For example, the mass and kinetic mixing will generate trilinear terms such as Z hZ, Z hA (with h being a CP-even Higgs boson) and Z W + W − , but these are suppressed by the small values of the Z −Z mixing, thus not changing the overall Z decay pattern significantly. In fact, we explicitly checked that these couplings change the overall bounds on the Z mass up to 5% at the most.
Having described the model and the relevant interactions, we can now present the collider bounds on the various Z realization we presented.

III. COLLIDER CONSTRAINTS
The main goal of this section is to derive the LEP and LHC bounds on the Z gauge bosons with charge assignments displayed in table I. We start by using LEP data [54] to constrain the kinetic and mass mixing terms which may significantly affect the Z properties if the mixing angle (ξ) is sufficiently large.
Typically this mixing angle is assumed to be arbitrarily small in models where there is no tree-level kinetic and mass mixing [55][56][57][58][59][60][61] but this assumption no longer applies to our models because we do include the kinetic mixing and the Higgs doublets are charged under the new gauge group (which implies the existence of mass mixing). In other words, once the gauge symmetry and the spontaneous EWSB mechanism are established, there is not much freedom left concerning the mass mixing, which is essentially set by the charges of the scalar fields under U (1) X and their VEVs. Moreover, the kinetic mixing, which can in principle be made ad hoc small, still arises at one-loop level since the SM fermions are charged under U (1) X [62][63][64]. After the analysis of the LEP constraints, we will then consider LHC bounds for generic Abelian extensions of the SM which, in the presence of sizeable kinetic mixing, have previously been obtained in [47,48,65].

A. LEP Limits
The LEP constraints on our models arise in the light of the excellent precision achieved by data collected and analyzed at such a machine. For example, by measuring processes such as e + e − → l + l−, where l = e, µ, LEP can restrictively probe beyond the SM scenarios that feature new particles coupling to charged leptons. In this connection, for our models, the presence of a massive Z that mixes with the Z leads to a deviation from the universal interactions of the latter with electrons and muons. This deviation can be parametrized as The parameter δΓ µe was measured at LEP [54] and can be presently used to place limits on the models studied in this work. Indeed the Z − Z kinetic and mass mixing will make this quantity depart from unit by inducing new interactions between the Z and SM fermions that are proportional to the mixing angle ξ, according to Eq. (8). Therefore, one can constrain such mixing using the LEP measurement of δΓ µe as shown in figure 1 by the dashed vertical lines. This idea was similarly exploited in the past to constrain different Z models [66][67][68][69][70][71][72][73].
We derived constraints on the mixing angle by plugging Eq. (8) into Eq. (10) and then comparing with the aforementioned LEP measurement [54]. In order to better understand these bounds, we recall that the mixing angle can be well approximated in the limit of large m Z by where G Xi are the couplings defined in Eq. (17) and encompass the gauge group dependence of the scalar doublets. Hence, using Eq. (11), we can relate the mixing angle ξ to the Z mass for all U (1) X models and thus find the lower mass bounds indicated by the vertical dashed lines in figure 1.
Once a gauge charge assignment is picked, G Xi is determined up to g X . Then, if we assume = 10 −3 and tan β = 10 (see Appendix A), in agreement with the latest experimental constraints on the 2HDM [74], the (sine of the) mixing angle sin ξ, g X and m Z are directly connected to one another leaving, in the end, only two free parameters. Their relation is model dependent and for this reason there is a curve for every model in figure 1. (The curve for model C is not visible in figure 1 because it is hidden between the curves for models D and B − L.) Notice that the vertical lines correspond to the lower mass bounds on the Z mass obtained by enforcing that all Z properties, crucially including its mass measurement m Z = 91.1876 ± 0.0021, remain in agreement with LEP data. We list these limits in table II as a function of the mixing angle ξ and Z mass.
It is now important to justify why our bounds do not explicitly depend on g X . Indeed, in figure 1 our limits rely on TeV   Table II. Summary of the LEP bounds on the (sine of the) mixing angle ξ vs m Z for all the U (1)X models in table I. the combination of sin ξ and m Z only. The dependence on g X enters in the Z mass and in sin ξ via G Xi , however, any change induced by g X can be parametrized as a shift on v S which is a free parameter and not an observable. In a nutshell, the Z mass can be approximated for large values of v S as where q X is the U (1) X charge of the SM singlet scalar. Then any change on the LEP limits due to g X can absorbed back by a redefinition of v S . From table II we conclude that Z masses below 2 TeV are excluded by LEP for all models. In particular, for the U (1) F model LEP imposes m Z > 4.25 TeV. We emphasize here that these limits result from Z − Z mixing effects and not from Z production at LEP. In Appendix A we also explicitly show how the Z mass changes with the mixing angle and which scale of symmetry breaking, v S , can be assumed in order to guarantee consistency with LEP data. Lower mass bounds of this nature are paramount for models that feature a large decay width lying outside the Narrow Width Approximation (NWA) which LHC limits are based on. We will come back to this point later on.
Anyhow, it is important to have these LEP bounds at hand since they already put strong limits on the mixing angles and on the Z masses for each of the models discussed in this work. In the derivation of such bounds we made some assumptions on the kinetic mixing and tan β. Nonetheless, the kinetic mixing is not much relevant because it only acts as a correction to the gauge coupling g X , see Eq.17, which, as already explained above, can be absorbed by a rescaling of v S . As for the value of tan β that enters in the Z mass, its impact is not important either because m Z is mainly set by Now that we have shown the LEP limits on our models we are ready to derive the LHC bounds.

B. LHC Limits
The best LHC bounds are obtained here by simulating at 13 TeV CM energy the process where f = q, l, leading to dijet or dilepton signals, respectively. Since this channel is mediated by a heavy Z , alongside γ and Z, a peak around the Z mass would appear at large values of the invariant mass of the dijet or dilepton final state. Initially, we will describe this in our simulation by adopting the aforementioned NWA, wherein the Breit-Wigner (BW) distribution capturing the propagation of the Z is replaced by a Dirac δ distribution. Eventually, we will allow for finite width effects as well. In the Monte Carlo (MC) generation we also need to account for the possible presence of up to two extra jets from QCD radiation alongside ff production, which are represented by X in the equation above. This obviously results into a better estimation of the production and decay cross section and its kinematics. What characterizes a dijet or dilepton signal in our study is of course the presence of an isolated jet or lepton pair with a large invariant mass. The numerical analysis is performed according to [39].
Signals with a resonant peak at high dijet or dilepton masses are of course absent in the SM model thus making the observation of such events a smoking gun signature for new physics, especially for a Z state. The main SM background stems from irreducible dijet and dilepton production via the γ and Z bosons, reducible tt production and decay as well as instrumental jet mis-reconstruction, but for invariant masses above 1 TeV they total a few events only for, e.g., a Center-of-Mass (CM) energy of 8 TeV and ∼ 20 fb −1 of integrated-luminosity [75].
In our work, firstly we aim at deriving LHC limits based on the dijet analysis performed by CMS with 13 TeV CM energy and L = 12 fb −1 [36] as well as the dilepton study conducted by ATLAS with L = 36 fb −1 [39], which are the corresponding most recent analyses for these datasets. To do so, we implemented all our U (1) X models in FeynRules [76] and simulated the partonic events with MadGraph5 [77]. We took into account hadronization and detector effects using Pythia8 [78] and Delphes [79], respectively, with the so-called k T -MLM jet matching scheme described in [80].
Since we are now discussing the on-shell production of a Z gauge boson (i.e., in NWA), the key quantities are: (i) the g X coupling that enters both in the production cross section and Z total width; the Z BR into (ii) light quarks and (iii) charged leptons; (iv) the Z mass; (v) the angle ξ that controls the Z − Z mixing. Notice that the latter is theoretically derived using Eq. (11) once g X and the Z mass are fixed. The fact that ξ is directly fixed by the model parameters makes our study more predictive. While its value is taken compliant with LEP data in all cases, we note that the LHC offers orthogonal and independent constraints on the models. We need to assess, however, which LHC experiment provides the most stringent bounds.
Regarding the Z BRs into light quarks and charged leptons, these change depending on the model under study because the fermions may have different quantum charges un-LHC bounds for gX = 0.1 Models Dilepton Table III. Dilepton bounds on our 2HDMs with U (1)X gauge symmetries at 13 TeV CM energy using L = 36 fb −1 for gX = 0.1. The dijet limits for gX = 0.1 are weaker in comparison thus are not displayed.
der U (1) X . The Z gauge boson features, in general, sizeable couplings to SM fermions making dijet and dilepton searches important environments to probe our models. Using the charge assignments defined above, we can therefore find lower bounds on the Z mass by comparing our predictions with the experimental limits on the overall production and decay rates.
The LHC exclusion plots in NWA in the aforementioned two channels are presented in figure 2. The CMS collaboration reported their limits in terms of the production cross section times BR into jets times acceptance (i.e., σ × BR × A, where A = 0.6 [36]) whereas ATLAS provides the bounds on the production cross section times BR into charged leptons only (i.e., σ × BR) [39]. We have also extrapolated the ATLAS experimental bound from 5 to 7 TeV making a leastsquares polynomial fit on the ATLAS data. Since at very large invariant masses errors are dominated by statistics, we can expect the experimental limit to indeed behave as reported on the right-hand side of figure 2, but we emphasize here that this extrapolation is not robust despite it provides a reasonable estimate of the ATLAS bound for Z masses above 5 TeV. Notice that, by comparing the two panels in figure 2, the dilepton limit is represented by a much smoother curve whereas the dijet one is rather bumpy due to the poorer reconstruction efficiency in the latter case. The first observation that we can easily draw by comparing the two plots in figure 2 is that dilepton bounds are more restrictive, as expected 1 .
These limits are displayed in table III for all models considered. We highlight that all these bounds are derived with g X = 0.1, i.e, of EW strength.
As we pointed out above, each Z couples differently to the SM fermions making then the total width Γ Z a modeldependent parameter. As Z couplings to fermions grow (i.e., g X increases), Γ Z also does, more so for large values of the Z mass, because of phase space effects. However, the relevant quantity playing a key role phenomenologically is the ratio between these two quantities, Γ Z /m Z . From now on, we will refer as NWA for the cases that accomplish Γ Z /m Z ∼ 1% and as Finite Width (FW) regime for the The solid red curve represents the observed 95% CL upper limit on the Z production cross section times BR times acceptance A for dijet events, where A = 0.6 (data are taken from Ref. [36]). Right panel: The solid red curve represents the observed 95% CL upper limit on the Z production cross section times BR for dilepton events (data are taken from Ref. [39]). The other curves account for the theoretical predictions as a function of m Z for all models discussed in this paper. The red dashed line in right panel represents the polynomial fit of the experimental data employed to extrapolate the ATLAS constraints up to 7 TeV. Here, gX = 0.1.
cases for which Γ Z /m Z ∼ 10%. It is important then to know when a particular model is in the NWA or FW setup.
In the Appendix we will show how Γ Z changes for each model while g X varies from 0.1 to 1, indeed covering both cases, NWA and FW, for all U (1) X models. Here, we have calculated in figure 3, for the BP with tan β = 10 and = 10 −3 , the single production cross section of Z times the dilepton BR for four cases of NWA and FW configurations, the ones without any cut on the dilepton invariant mass (NWA and FW) and the ones with the "magic cut" of [81], i.e., |m ll − m Z | < 5% √ s (NWA-MC and FW-MC), with √ s = 13 TeV. Recall that Ref. [81] adopted such a constrain to capture interference effects between Z signals and corresponding SM irreducible backgrounds in such a way that these are minimized over the relevant kinematic range, thus enabling one to perform quasi-model-independent analyses, essentially preserving the NWA scheme most often used by the experimental collaborations for Z boson searches. Let us stress that each U (1) X model attains NWA or FW approximations for different values of g X , for example, the U (1) G model achieve FW conditions for g X = 0.8 while the U (1) F model does so for g X = 0.3.
The first observation that we can draw from figure 3 is that the magic cut has no impact at all on the cases in NWA, as expected, and a rather moderate impact in the FW regime and only for large masses, where interference effects are more substantial, the effects being essentially the same for all models considered. In this respect, it is worth noting that, since the FW regime is achieved here by increasing the coupling g X , this implies not only a larger width but also a much larger cross section. In fact, in this case, at large dilepton invariant mass, the resonant Z contribution is dominant over the interference between the Z contribution and the SM one due to γ and Z exchange, as generally g X > e and e/ sin θ W plus the Z is resonant while γ and Z are highly off-shell. In these conditions, the NWA result and the one in FW are very similar, as already shown in Ref. [82]. We can then use the ATLAS data of Ref. [39] which assume a narrow Z , even when Γ Z /m Z = 10%, which is here obtained by suitably adjusting g X (differently) for all models considered in order to extract limits. We do so in the left frame of figure 4. We observe that in this FW case the excluded m Z range is comparable to the one of the previous case in NWA, as expected.
The magic cut becomes important when we obtain the prescribed Γ Z /m Z = 10% ratio by enforcing the latter by hand while adopting g X = 0.1 throughout our models. This is well justified. Recall in fact that U (1) X models like those considered here can be remnants of Grand Unification Theories (GUTs), wherein there can be large particle spectra (of both fermions and bosons) into which such heavy Z s can decay to, so as to justify a 10% width-to-mass ratio. Besides, in experimental searches, mass and width are treated as independent parameters. In this case, interference effects are substantial so that the constraint advocated by Ref. [81] is mandatory. We thus need to use experimental data selected using the magic cut, which are those of Ref. [83]. In this analysis, see figure  4 therein, CMS plots the upper limits at 95% CL on the product of production cross section and BR for a Z with finite width, including the case of 10% of the resonance mass, relative to the product of production cross section and BR for the Z boson. Hence, we present the right frame of figure 4. We observe here that the discussed interference effect is substantial, as a much more restricted m Z range can be probed with respect to the cases when this is negligible, signalling that the correction is predominantly negative.
We can now compare our LHC bounds with those from LEP and notice that they are rather complementary to each other. For instance, LHC requires m Z > 3.6 TeV for g X = 0.1 for the U (1) F model even when away from the NWA, see Appendix, though this is subject to systematic errors (of modeldependent FW and interference effects that ought to be quantified for each energy and luminosity values). LEP instead  Figure 3. Production cross section of a Z times its dilepton BR for each model considered assuming a NWA (Γ Z /m Z ∼ 1%) and FW regime (Γ Z /m Z ∼ %10). Alongside the rates without any cuts we also display those following the magic cut of [81], i.e., |m ll − m Z | < 5% √ s, which are denoted by the NWA-MC and FW-MC labels. [σ Figure 4. Left panel: The solid red curve represents the observed 95% CL upper limit on the Z production cross section times BR for dilepton events (data are taken from Ref. [39]) while the other curves are for the theoretical predictions as a function of m Z taking different gX values yielding Γ Z /m Z = 10% for all models considered. Right panel: The solid red curve represents the observed 95% CL upper limit on the Z production cross section times BR relative to the Z production cross section times BR for dilepton events (data are taken from Ref. [83]) while the other curves are for the theoretical predictions as a function of m Z taking gX = 0.1 and setting by hand Γ Z /m Z = 10% for all models considered. The red dashed line in both panels represents the polynomial fit of the experimental data employed to extrapolate the ATLAS and CMS constraints up to 7 TeV. Notice that the magic cut of [81], i.e., |m ll − m Z | < 5% √ s has been applied here on both CMS data and MC predictions (right panel only).
sets a lower mass limit of m Z > 2.5 TeV constituting a weaker but more robust bound on the Z mass on this case. This clearly demonstrates that LEP data is still a powerful tool to constrain new physics models, especially those that feature a large width-to-mass ratio. Such scenarios are commonly used in dark matter model building endeavors [84], for example. For model U (1) D , in contrast, the Z does have a narrow width for g X = 0.1, with LHC(LEP) excluding Z masses below 3.5 TeV(2.8 TeV). In this case, the simpler NWA analysis described above is applicable and the LHC searches provide not only the strongest bound on the Z mass but also a very solid one.
In summary, one needs to truly consider both LEP and LHC data to extract reliable bounds on the Z mass and couplings of U (1) X models, bearing in mind that these limits are complementary and orthogonal to each other.

IV. HL-LHC AND HE-LHC SENSITIVITY
Bearing in mind that the models studied in this work predict a large Z width, relative to its mass, for g X ∼ 1 and that this may imply a reassessment of the validity of the magic cut in presence of much larger luminosities and/or energies of the LHC, including that of systematic effects from both the theoretical and experimental side, we will derive the projected sensitivity for the HL-LHC and HE-LHC configurations only for g X = 0.1.
The HL-LHC setup is characterized by L = 300 and 3000 fb −1 while the HE-LHC represents the LHC upgrade phase with CM energy of 27 TeV. To find their physics sensitivity we will adopt the strategy described in [85]. In short, what the code does is to solve an equation for M new (i.e., the new limit on m Z ), knowing the current bound M , as follows: with obvious meaning of the subscripts. The results from this iteration are summarized in table IV. The lower mass bounds found presently compared to the expected ones at the HL-LHC and/or HE-LHC clearly show how important is any LHC upgrade to test new physics models including a Z . For some models such as, e.g., the U (1) A , the HE-LHC will potentially probe Z masses up to 7 TeV, i.e., three times higher than with current data. In summary, the HL-LHC and (especially) the HE-LHC will represent discovery machines being able to probe gauge boson masses with unprecedented sensitivity up to the O(10 TeV) domain.

V. CONCLUSION
2HDMs offer an interesting framework for Higgs boson phenomenology above and beyond what offered by the SM but, at the same time, also face the presence of dangerous FC-NCs which are severely constrained by data. An ad-hoc Z 2 symmetry is usually imposed to prevent one of the two scalar doublets from generating fermions masses and thus freeing 2HDMs from these constraints. In this work, we derived collider bounds on 2HDMs that instead address the flavor problem via gauge symmetries and naturally accommodate neutrino masses via a type-I seesaw mechanism. These gauge symmetries give rise to Z gauge bosons that feature different coupling structures to leptons and quarks. We thus exploited the complementarity between dijet and dilepton data from the LHC to find lower mass bounds on the corresponding Z gauge bosons, by investigating scenarios where these are applicable in both NWA and FW regime, finally contrasting the ensuing limits extracted from LEP data which are sensitive, on the other hand, to the effects of the kinetic and mass Z − Z mixings. Lastly, we have presented the sensitivities of the HL-LHC and HE-LHC to such new physics, showing that they are capable of probing Z masses up to a factor three higher than presently. .

ACKNOWLEDGMENTS
The authors thank Werner Rodejohann and Pyung-won Ko for discussions and comments. DC and FSQ acknowledge financial support from MEC and UFRN. FSQ also acknowledges the ICTP-SAIFR FAPESP grant 2016/01343-7 for additional financial support. SM is supported in part by the NExT Institute and acknowledges partial financial support from the STFC Consolidated Grant ST/L000296/1 and the H2020-MSCA-RISE-2014 grant no. 645722 (NonMinimal-Higgs).

APPENDIX A
We collate here some information aiding the understanding of the main part of the paper, by dividing it into sections relating to key computational aspects.

GAUGE BOSON COUPLINGS
After rotating to a basis in which the gauge bosons have canonical kinetic terms, the covariant derivative in terms of small reads or, explicitly, where we defined for simplicity with Q Yi being the hypercharge of the scalar doublet, which in the 2HDM is taken equal to +1 for both scalar doublets, and Q Xi is the charge of the scalar doublet i under U (1) X .
Then the part of the Lagrangian responsible for the gauge boson masses becomes where v 2 = v 2 1 + v 2 2 . Eq. (18) can then be written as with and Summarizing, after the symmetry breaking, one can realize that there is a remaining mixing between Z 0 µ and X µ that can be expressed through the symmetric matrix or, explicitly, The above expression, Eq. (24), representing the mixing between the Z 0 µ and X µ bosons, is given as function of arbitrary U (1) X charges of doublet (or singlet) scalars. It is important to notice that, when Q X1 = Q X2 and there is no singlet contribution, the determinant of the matrix Eq. (24) is zero.
The matrix in Eq. (24) is diagonalized through a rotation O(ξ) and its eigenvalues are The ξ angle is given by Since this mixing angle it supposed to be small, as m 2 Z m 2 Z , we can use tan ξ ∼ sin ξ with We can expand this equation further to find a more useful expression. Substituting the expressions for G Xi and factoring out the m Z mass, we finally get sin ξ m 2 Z m 2 Z g X g Z (Q X1 cos 2 β + Q X2 sin 2 β) + tan θ W .

(29)
Z PROPERTIES AND LEP BOUNDS All the U (1) X configurations demand high v S > ∼ 5 TeV values in order to be in agreement with the SM Z mass measurements of m Z = 91.1876 ± 0.0021. This has been taken into account and the LEP limits showed here encompass this constraint. We compiled them in table II to ease the reading. In this procedure we adopted = 10 −3 to be consistent with EW precision data [86] and tan β = v 2 /v 1 = 10 in agreement with the latest experimental constraints on the 2HDM [74]. In the left panel of figure 5 we showed how the Z mass changes depending on the mixing angle adopted for every single U (1) X model. This illustrates that, indeed, only small mixing angles of the order of 10 −3 are allowed by LEP data. In the right panel of the same figure, for completeness, we exhibit a heat map to indicate which value of v S is needed to reproduce a small mixing. Using LEP precision data and the theoretical predictions discussed above, we can estimate that v S should be larger than ∼ 10 TeV for all models investigated here.

Z WIDTH
In figure 6 we show how the Z mass scales with v S for g X = 0.1 and 1. It is clear that, irrespectively of the details of the model, as soon as v 2 S 246 2 GeV, the Z mass is dominated by v S and scales linearly with it. In figure 7 we show the Z width as a function of its mass for g X = 0.1, 0.4, 0.7 and 1, where one can easily see that, for increasing g X , the Z width eventually becomes too large relative to the mass so that one cannot use the NWA to derive LHC bounds on the Z properties.