Loopholes in $Z^\prime$ searches at the LHC: exploring supersymmetric and leptophobic scenarios

Searching for heavy vector bosons $Z^\prime$, predicted in models inspired by Grand Unification Theories, is among the challenging objectives of the LHC. The ATLAS and CMS collaborations have looked for $Z^\prime$ bosons assuming that they can decay only into Standard Model channels, and have set exclusion limits by investigating dilepton, dijet and to a smaller extent top-antitop final states. In this work we explore possible loopholes in these $Z^\prime$ searches by studying supersymmetric as well as leptophobic scenarios. We demonstrate the existence of realizations in which the $Z^\prime$ boson automatically evades the typical bounds derived from the analyses of the Drell-Yan invariant-mass spectrum. Dileptonic final states can in contrast only originate from supersymmetric $Z^\prime$ decays and are thus accompanied by additional effects. This feature is analyzed in the context of judiciously chosen benchmark configurations, for which visible signals could be expected in future LHC data with a $4\sigma-7\sigma$ significance. Our results should hence motivate an extension of the current $Z^\prime$ search program to account for supersymmetric and leptophobic models.


Introduction
Although the discovery of the Higgs boson at the LHC has allowed for the completion of the particle spectrum of the Standard Model (SM), the issue of its extension still stands. Despite the experimental success in predicting most data observed so far, the SM indeed exhibits several limitations and shortcomings that motivate the study of beyond the Standard Model theories. Among those, supersymmetry, and in particular its minimal incarnation known as the Minimal Supersymmetric Standard Model (MSSM), is one of the most appealing options. Unifying internal and external symmetries, supersymmetry provides a natural solution to the long-standing hierarchy problem, allows for gauge-coupling unification at high energies and predicts a stable particle that could address the problematics of Dark Matter. Despite these numerous motivations, no compelling evidence for supersymmetry has been found and the MSSM starts to be heavily constrained. Moreover, the MSSM suffers from severe fine-tuning issues related to the discovery of a SM-like Higgs boson, as well as the lack of any satisfactory explanation for the magnitude of the supersymmetric bilinear Higgs mass parameter µ that must unnaturally be of the order of the electroweak symmetry breaking scale.
As a consequence, arguments have been raised in favor of extending the MSSM superfield content by at least one singlet chiral superfield. Its scalar component can induce both supersymmetry breaking and dynamical generation of the µ term by getting a non-vanishing vacuum expectation value at the minimum of the scalar potential [1][2][3][4][5]. Such singlet superfields also appear under supersymmetric scenarios where the Standard Model gauge group is extended, the scalar singlet yielding the breaking of the additional gauge symmetry [6,7]. This setup is furthermore motivated in a grand-unified scheme where a restricted set of high-dimensional representations are used to encompass all MSSM supermultiplets, and where all gauge couplings unify. In this context, the necessity of using representations of the unified gauge group automatically leads to the introduction of right-handed neutrino superfields, which consequently provides a solution for neutrino-mass generation, as well as vector-like fermions.
Among all Grand Unification Theories (GUT), those based on gauge groups of rank 6, named E 6 , have been extensively discussed as interesting possibilities [7,9,10]. In particular, the breaking pattern of E 6 to the electroweak symmetry results in the appearance of extra U (1) symmetries. From a bottom-up perspective, extending the MSSM with the introduction of an extra U (1) gauge group has numerous advantages, namely forbidding a too rapid proton decay without introducing an ad hoc discrete R-parity symmetry and making all field masses stable with respect to quantum corrections. Moreover, it is always possible to choose the U (1) field charges to ensure anomaly cancellation and gauge-coupling unification. Besides, the U (1) models do not suffer from the presence of cosmological domain walls, unlike theories like the Next-to-Minimal Supersymmetric Standard Model [11,12].
While the spectrum of U (1) supersymmetric (UMSSM) models is altered from that of the MSSM, the most secure prediction emerging from the extended gauge symmetry consists of the existence of a novel neutral Z boson, like in non-supersymmetric U (1) extensions of the SM. This makes the Z boson a prime target for the LHC physics program, as the proof of its existence would constitute a promising indicator of a more general gauge structure. Any gauge group of rank greater than four (any group larger than SU (5)) indeed leads to the appearance of at least one extra neutral gauge boson. All current Z analyses at the LHC are however guided by non-supersymmetric considerations in which the Z boson only decays into SM particles [13][14][15][16][17]. Besides E 6 -inspired Z , the experimental collaborations have also explored the so-called Sequential Standard Model (SSM), the simplest extension of the Standard Model, wherein Z and possible W bosons have the same couplings to fermions as the Z and W . This model is not theoretically motivated, but it is often used as a benchmark for the analyses, since the production cross section in the SSM just depends on the extra boson masses.
Along these lines, the ATLAS and CMS collaborations have searched for Z bosons by investigating dilepton and dijet final states. In detail, by using high-mass dilepton data at 13 TeV, the ATLAS collaboration [18] set the mass exclusion limits M Z > 4.5 TeV in the SSM and M Z > 3.8-4.1 TeV in U(1) models, whereas CMS obtained M Z > 4.0 TeV (SSM) and M Z > 3.5 TeV (GUT-inspired models) [19]. For dijets, the limits are much milder and read M Z > 2.1-2.9 TeV (ATLAS) [20] and M Z > 2.7 TeV (CMS) [21].
In a UMSSM framework, the inclusion of the supersymmetric decay modes of the Z bosons may nonetheless change these conclusions [22][23][24][25][26][27]. Above all, the opening of new decay channels lowers the branching ratios into SM final states and therefore the Z mass exclusion limits. In fact, Ref. [28] found an impact of about 200 GeV on the mass exclusion limits by comparing the 8 TeV ATLAS and CMS data on high-mass dileptons with UMSSM predictions for a benchmark point of the parameter space. Furthermore, in the UMSSM, a leptophobic Z can yield the production of dilepton final states only through cascade decays into intermediate electroweakinos, which contrasts with the leptophobic non-supersymmetric case where this is simply not allowed [29]. The bounds on the Zboson mass and production cross section derived from the above-mentioned searches should then be revisited when more general theoretical contexts like the UMSSM or leptophobia are considered.
On different grounds, the hadronic environment at the LHC is so complex that new physics searches always rely on some simplifying assumptions in the form of the potential signals. For instance, most supersymmetry searches have been designed from the idea on how the MSSM could manifest itself in a typical LHC detector: they may hence be not suitable for given non-minimal supersymmetric realizations. In the UMSSM framework, which we focus on in this work, we consider Z -boson signals that can potentially differ from the non-supersymmetric case. We restrict our analysis to leptonic Z decay modes that are easier to explore, even if the expected signals are plagued by larger SM backgrounds. We additionally focus on UMSSM realizations in which the Z boson is leptophobic, but where it could give rise to leptonic signatures via supersymmetric cascade decays into leptons and missing energy. This therefore offers an alternative opportunity to find both an extra gauge boson and supersymmetry from the study of the decays of a resonantly-produced colorless particle. This is one of the scenarios that we wish to investigate in this work, after imposing the most up-to-date constraints on the model. We hence aim at providing a clear roadmap for the discovery of unconventional leptophobic Z bosons, such as those that could arise in UMSSM scenarios and that escape detection when only considering standard LHC searches for extra gauge bosons.
Our work is organized as follows. In Section 2, we briefly introduce U (1) supersymmetric models as when the gauge symmetry is designed as emerging from the breaking of an extended E 6 symmetry at the grand unification scale. We pay particular attention to the mass, mixing patterns and interactions of the extra neutral gauge boson and show under what conditions it could be made leptophobic. We finally set up the parameter-space region to be scanned over and proceed to its exploration in Section 3, focusing on two different way to impose boundary conditions. In Section 4, we concentrate on scenarios where the Z boson does not directly decay into leptons and study its phenomenology at colliders, highlighting a preferred selection strategy that could lead to its discovery. We summarize our results and conclude in Section 5.

Theoretical framework
There are different ways to implement a U (1) extension in the MSSM: one of the most commonly used parameterizations is inspired by grand-unified models, based on a rank-6 group E 6 , where the symmetry-breaking scheme proceeds via multiple steps,  The U (1) symmetry that survives at the electroweak scale is taken as a linear combination of U (1) χ and U (1) ψ , where we have introduced the E 6 mixing angle θ E 6 . The neutral vector bosons associated with the U (1) ψ and U (1) χ gauge groups are called the Z ψ and Z χ bosons, while a generic Z is given by the mixing of these Z ψ and Z χ states, as in Eq. (2.2). Different U (1) models can be classified according to the sole value of the θ E 6 mixing angle, and the charges Q of the supermultiplets are fixed to ensure the theory to be anomaly-free. Six popular setups are summarized in Table 1, with the corresponding Q charges listed in Table 2. In the notations of this last table, Q and L denote the left-handed weak doublets of quark and lepton fields, H u and H d the two weak doublets of Higgs fields, U and D the right-handed weak singlets of up-type and down-type quarks, E and N the right-handed weak singlets of charged leptons and neutrinos, and S a scalar singlet. In the case of supersymmetric extensions of the Standard Model, such as the MSSM, all fields in Table 2 must actually be understood as superfields containing also the supersymmetric partners of the fermions and Higgs bosons. In principle, the matter sector of E 6 should also feature vector-like exotic (s)quarks Q D andQ D which have the same U (1) charges as the H u and H d fields, respectively [9]. In the following, we assume that these exotic states are too heavy to be relevant at LHC energies and neglect them in our phenomenological analysis 1 .
The Higgs supermultiplet content (H u , H d and S) is large enough to allow both for the breaking of U (1) via the scalar singlet field s, and of the electroweak symmetry through the neutral components of the scalar Higgs doublets h u and h d . All electrically-neutral Higgs fields indeed get non-vanishing vacuum expectation values at the minimum of the potential and carry non-trivial U (1) charges.
In the grand-unified framework, the field content is organized into vector representations (27) of the E 6 group; the latter further branches as 27 = 16 ⊕ 10 ⊕ 1 into the irreducible representations of the SO (10) subgroup that arises at the first step of the E 6 breaking scheme of Eq. (2.1). In the conventional field assignment, the representation 16 contains the left-handed quark and lepton supermultiplets (Q and L), as well as the righthanded quarks and leptons (U , D, E and N ), while the Higgs fields (H u and H d ) and the exotic quarks Q D andQ D are in the representation 10. An alternative framework consists of having instead H u andQ D lying in the 16 and L and D in the 10 representation. According to whether one chooses the standard or unconventional assignment, the phenomenology of the Z boson may be different. In the following, we shall adopt the standard SO(10) representation choices, with the exotic quarks lying in the 10 representation. Nevertheless, the unconventional scenario can be easily recovered by redefining θ E 6 → θ E 6 + arctan √ 15 in Eq. (2.2) [30].
In principle, the Higgs fields in the 27 representation of E 6 should occur in three generations. However, as discussed in Refs. [8,31,32], it is always possible to perform a unitary transformation to a basis where only one generation of Higgs bosons gets a nonvanishing vacuum expectation value. The scalars with zero vacuum expectation values were called 'unHiggs' in Refs. [8,31]. Through our analysis, we shall neglect the two generations of such states and focus on the 'true' Higgs bosons, which exhibit a non-zero vacuum expectation values and are denoted by H u and H d .
The 16 representation of SO(10) is then decomposed in terms of those of SU (5) as 16 = 10⊕5⊕1. The 10 representation of SU (5) is suitable to include right-handed up-type quark and charged-lepton supermultiplets, together with the weak doublets of left-handed quarks, whereas the 5 representation contains right-handed down quarks and left-handed lepton supermultiplets; the 1 representation includes right-handed (s)neutrinos [33]. The UMSSM superpotential is thus given, all flavor indices being omitted for clarity, by: The Yukawa interactions are encoded in a set of four 3 × 3 matrices in flavor space, Y u , Y d , Y l and Y ν , and the strength of the supersymmetric Higgs self-interactions is described by the λ parameter. After the breaking of the U (1) symmetry, this λ-term induces the dynamical generation of an effective µ-term (denoted µ eff in the following) that allows for the resolution of the so-called MSSM µ-problem 2 . As supersymmetry has to be softly weak isospin quantum numbers allowing for a superpotential interaction term involving ordinary quarks and inducing rapid proton decay. Their mass must thus be comparable to the GUT scale to prevent the proton from decaying too quickly [9].
2 µ eff is related to λ and to the vacuum expectation value of the scalar singlet s via µ eff = λ s . broken, we introduce in the Lagrangian explicit mass terms for all gaugino and scalar fields, where the U (1) Y , U (1) , SU (2) L and SU (3) c gaugino Weyl fermions are denoted by λB, λB , λW and λg, respectively, and where h d , h u , s,q,d † ,ũ † ,l,ẽ † andν † are the scalar components of the H d , H u , S, Q, D, U , L, E and N superfields. The set of M i and m i parameters moreover denote the soft gaugino and scalar mass parameters, respectively. Additional soft terms, related to trilinear scalar interactions, are also present and can be derived from the structure of the superpotential, where the A e , A ν , A d and A u 3×3 matrices stand for the strengths of the soft Higgs-boson interactions with charged sleptons, sneutrinos, down-type squarks and up-type squarks, respectively. The A λ parameter is finally related to the trilinear soft multi Higgs-boson coupling.
In order to calculate the sfermion masses, one would need to set up an explicit framework for supersymmetry breaking, such as a gauge-, gravity-or anomaly-mediated mechanisms, which goes beyond the goals of the present paper. We only recall that supersymmetry can be spontaneously broken if the so-called D-term and/or F -term in the scalar potential have non-zero vacuum expectation values. The F -terms are proportional to the SM particle masses, and are therefore important only for stop quarks, whereas D-terms are relevant for both light and heavy sfermions and contain contributions due to electroweak symmetry breaking and, in case of extension of the MSSM, to the Higgs bosons which break the extended symmetry [22,25,26]. Hereafter, we account for F -and D-term corrections to the sfermion masses, but do not present their explicit expressions, for the sake of brevity.
After the spontaneous breaking of the symmetry group down to electromagnetism, the W , Z and Z bosons get massive and the photon stays massless. In general, for a U (1) extension of the SM, there is mixing between the Z and Z eigenstates, parameterized by a mixing angle α ZZ . However, electroweak precision data strongly constrain α ZZ to be very small [34]. At tree level, the squared masses of the Z and Z bosons are given by: where h 0 d and h 0 u stand for the neutral components of the down-type and up-type Higgs fields h d and h u and g 1 , g 2 and g are the coupling constants of the U (1) Y , SU (2) L and U (1) gauge groups, respectively. As discussed, e.g., in Ref. [9], whenever the singlet s has a large vacuum expectation value (which contributes only to the Z mass), as will be the case hereafter, M 2 In the Higgs sector, as discussed above, one should deal with three generations of Higgs fields, although, in our chosen basis, only one generation (the so-called 'true' Higgs bosons) exhibits non-zero vacuum expectation values. Mass mixing matrices and mass eigenstates of the two generations of Higgs bosons with zero vacuum expectation values are thoroughly debated in [31]. In principle, because of the presence of these other states, one should impose further constraints on our scenario coming, e.g., from the current measurements of the (SM-like) neutral-Higgs production cross section and branching ratios, as well as from the exclusion limits on charged-Higgs bosons. In our work, however, such extra Higgs states and related constraints will be neglected.
In fact, after electroweak symmetry breaking, for each generation of Higgs fields, one is left with two charged and four neutral scalar bosons, namely one pseudoscalar and three neutral scalars, including a novel singlet-like scalar Higgs, inherited by the U (1) symmetry. In the following, we shall account for only one generation of Higgs bosons and denote by H ± the charged bosons, h and H the MSSM-like neutral scalars, with h roughly corresponding to the Standard Model Higgs, A the pseudoscalar and H the extra scalar associated with the U (1) gauge group.
As discussed, e.g., in Ref. [8], for s much larger than h 0 u and h 0 d , diagonalizing the neutral Higgs mass matrix is straightforward and the singlet-like H has mass M 2 H g 2 Q S 2 s 2 , hence it is roughly degenerate with the Z , according to Eq. (2.6). The other neutral Higgs H has instead approximately the same mass as the pseudoscalar A and as the charged H ± : as a result, the heaviest scalar Higgs of the spectrum could be either H or H , depending on whether the Z is lighter or heavier than A.
In the gaugino sector, with respect to the MSSM, one has two extra neutralinos, related to the supersymmetric partners of Z and H bosons, which yields a total of sixχ 0 1 , . . . ,χ 0 6 neutralino states. As discussed in Ref. [26], the newχ 0 5 andχ 0 6 eigenstates are often too heavy to contribute to the Z phenomenology at the LHC. As the new Z is electrically neutral, the chargino sector stays instead unchanged with respect to the MSSM.
On top of mass mixings, both U (1) Y and U (1) bosons are allowed to mix kinetically [35]. The corresponding Lagrangian reads, in terms of the gauge boson component fields, whereB µν andẐ µν are the U (1) Y and U (1) boson field strength tensors, respectively, and χ is the kinetic mixing angle. In order to understand the physical implications of the kinetic mixing, it is necessary to diagonalize the field strengths, which is achieved via a GL(2, R) rotation, whereB µ and Z µ are the original U (1) and U (1) gauge fields, with non-diagonal kinetic terms, while B µ and Z µ have now canonical diagonal kinetic terms. As discussed in Refs. [9,35], for M 2 Z M 2 Z and small values of χ, the impact of the kinetic mixing on the gauge boson masses is negligible. It nonetheless can have a significant effect on the coupling of the Z boson with fermions. In fact, the interaction Lagrangian of the fieldsB µ andẐ µ with a generic fermion ψ i , with charges Y i and Q i under the U (1) and U (1) groups, is given by which can then be rewritten in terms of B µ and Z µ as Leptophobic scenarios can hence be obtained requiringQ L =Q E = 0 [36][37][38]. Since Y L = −1/2 and Y E = 1, Eq. (2.11) dictates that leptophobia can be achieved only if Q E = −2Q L : this relation between the doublet and singlet leptonic charges is typical for the U (1) η configuration, as shown in Table 2. Furthermore, if one assumes, as will be done in the following, the typical GUT-inspired relation between the U (1) and U (1) couplings g 1 (M Z )/g (M Z ) = 3/5, then leptophobia requires the additional condition sin χ ≈ −0.3.
As a result, we expect leptophobic Z models to naturally arise for E 6 mixing angles in the neighbourhood of θ E 6 θ η ± nπ, n = 0, 1, 2, 3, . . . , (2.12) with the Z -boson leptonic couplings being either exactly zero or very suppressed. In the following, we shall account for the kinetic mixing of U (1) Y and U (1) gauge groups, with the U (1) charges of all our matter fields given by Eq. (2.11).

Parameter-space scan and constraints
UMSSM theories rely on numerous free parameters so that simplifying assumptions are in order for a practical parameter-space exploration. Hereafter, we impose minimal flavor violation, so that all the flavor-violating parameters of the soft supersymmetry-breaking Lagrangian are considered as vanishing, and enforce unification boundary conditions on the remaining soft parameters.
In the first class of scenarios which we investigate, unification is assumed to occur at a very high scale M GUT ≈ O(10 16 ) GeV and all parameters are then run down to M Z according to renormalization group evolution. More precisely, all gauge couplings are assumed to unify at a given high scale and the U (1) coupling is enforced to satisfy Furthermore, all scalar masses are set to a common value M 0 , whilst all gaugino masses are taken equal to another universal mass M 1/2 . All trilinear soft couplings are assumed to be proportional to the respective Yukawa coupling matrices with a universal proportionality factor A 0 , so that (2.14) Parameter Scanned range Parameter Scanned range Parameter Scanned range Parameter Scanned range In the Higgs sector, we fix the values of the effective µ eff parameter, the ratio of the vacuum expectation values of the neutral components of the two Higgs doublets tan β = v u /v d , the trilinear soft coupling A λ , as well as the Z mass M Z . Finally, the diagonal entries of the neutrino Yukawa coupling matrices are set to a very small value, O(10 −11 ), in such a way as to ignore the sneutrino soft trilinear interactions. The ensemble of free parameters considered in our exploration of the UMSSM parameter space is thus given by where we have additionally included the E 6 mixing angle θ E 6 . We vary those parameters over the ranges given in the top panel of Table 3.
In the second class of scenarios considered in this work, unification is imposed at the Z mass scale. In this case, we just enforce the unification of the trilinear couplings as in Eq. (2.14) and set all scalar and gaugino masses being kept free. The entire set of free parameters is thus here given by with the ranges over which those parameters vary presented in Table 3. In our scanning procedure, we analyze all possible anomaly-free UMSSM models derived from the breaking of an E 6 gauge symmetry. We generate the particle spectrum by making use of the Sarah code, version 4.6.0 [47], and its interface to SPheno 3.3.8 [48]. In order to test the phenomenological viability of the model, we compute various properties of the Higgs sector, such as the mass of the lightest Higgs state and the corresponding collider signal strengths by means of the HiggsBounds (version 4.3.1) and HiggsSignals Mτ > 81 GeV [42] Mẽ > 107 GeV [42] Mμ > 94 GeV [42] Mt > 900 GeV [43] BR Table 4. Experimental constraints imposed within our scanning procedure in order to determine the parameter-space regions of interest.
(version 1.4.0) packages [49,50]. The scan itself and the numerical analysis performed in this work have been achieved by interfacing all programs using also the pySLHA package, version 3.1.1 [51]. The parameter space is probed by using the Metropolis-Hasting sampling method, requiring consistency with the experimental bounds on masses and decay rates shown in Table 4. In particular, we require the mass of the Standard Model Higgs boson to agree with the measurements up to an uncertainty of 3 GeV, and the χ 2 fit of the available Higgs signal strengths is bounded to be smaller than the conservative value of 70. Other constraints, connected to the bounds on the masses of supersymmetric particles and on several flavor observables, are evaluated relying on the SPheno code. This includes in particular tests of the strict limits stemming from B-meson decays [43][44][45]. As for the supersymmetric sector, we enforce the LEP limits on slepton, chargino, and neutralino masses quoted in Ref. [42], while for gluinos and stops we implement the bounds set by CMS [41] and ATLAS [43], respectively.

Supersymmetric Z Phenomenology
In this section, we analyze the phenomenology of the two classes of UMSSM scenarios introduced in Section 2.2. In the subsequent Section 4, specific configurations where the Z boson is leptophobic by virtue of the kinetic mixing of U (1) Y and U (1) are in contrast investigated.
In order to apply the LHC constraints on the properties of Z bosons, we calculate the Z production cross section at next-to-leading order (NLO) accuracy in QCD [52,53]. This relies on the joint use of FeynRules version 2.3.27 [54] and the included NLOCT package [55], as well as FeynArts [56], for the automatic generation of a UFO library [57] containing both tree-level and counterterm vertices necessary at NLO. This UFO model is then used by MadGraph5 aMC@NLO (version 2.5.5) [58] for the numerical evaluation of  Table 5. g values and dilepton branching ratios for commonly studied U (1) models with UMSSM parameters satisfying the constraints detailed in subsection 2.2. Quoted are g min , the minimum value of g (M Z ), along with the corresponding spread ∆g and the smallest possible branching ratio into leptons with (UMSSM) and without (USM) supersymmetric contributions to the Z decays. the hard-scattering matrix elements, which are convoluted with the NLO set of NNPDF 2.3 parton distribution functions (PDF) [59]. Using the decay table provided by the SPheno package and assuming the narrow-width approximation, we compare our predictions with the ATLAS limits on Z bosons in the dilepton mode [18] in order to estimate the impact of supersymmetric decay channels.

Scenarios With High-Scale Boundary Conditions
In this subsection, we focus on our first class of UMSSM scenarios where the proportionality between g and g 1 is imposed at the GUT scale and where all free parameters in Eq. (2.15) are fixed at M GUT and then evolved down to the Z scale by means of renormalization group equations.
We have found that some parameter regions satisfying the constraints in Table 4 exist for a wide set of values of the E 6 mixing angle θ E 6 . The LHC collaborations typically use the rate σB ≡ σ(pp → Z ) × BR(Z → l + l − ) to obtain the exclusion limits on the Z mass. For the sake of exploring possible loopholes in the Z searches, we are therefore especially interested in scenarios which minimize the σB product, namely featuring small values of the g coupling and of the BR(Z → l + l − ) branching ratio. In fact, when running the renormalization group equations, scanning the parameters in the ranges presented in Table 3, imposing the constraints of Table 4 and accounting for proper threshold matching conditions, g (M Z ) ends up with lying in a range [g min , g max ].
In Table 5 we quote, for a few U (1) models, the minimum value of g at the M Z scale and the spread ∆g , defined as and expressed as a percentage. The minimum branching fraction of Z decays into dilepton final states, including supersymmetric channels (UMSSM) and without supersymmetry (USM) is also quoted. In the table, we have discarded the models U (1) χ and U (1) S . As discussed, e.g., in Refs. [25,27], U (1) χ models are ill-defined in supersymmetry as it they typically lead to unphysical sfermion masses after adding to the soft masses the D-term contributions. As to U (1) S , it may be theoretically acceptable, but we were not able to find scenarios capable of satisfying the constraints of Table 4. From Table 5, we learn that the deviations of g from g min are rather small, with ∆g being of at most about 8%, but the impact of the inclusion of supersymmetric decays on the dilepton branching fraction is remarkable for most models. In the U (1) ψ and U (1) η scenarios, for example, BR(Z → ll) decreases by about 35% and 25%, respectively, once decays into sfermions and gauginos are accounted for. Nevertheless, all models still exhibit substantial dilepton Z decay rates, varying between 3% and 10%.
In Fig. 1 we compare the ATLAS limits on high-mass dileptons at the 1σ (green) and 2σ (yellow) levels with our predictions for σB, obtained in the context of U (1) ψ and U (1) η (upper panel), as well as U (1) I and U (1) N (lower panel) gauge groups, in the range 2 TeV < M Z < 5 TeV. We consider both supersymmetric (markers with error bars) and non-supersymmetric cases (dashed lines) and include NLO QCD corrections to the production cross section σ(pp → Z ). The error bars around the supersymmetric results include two contributions: first, they account for the spread covered in the scan and second, they include the theoretical uncertainties stemming from traditional scale and parton density variations in the NLO computation. We found that the latter uncertainty varies from 5% for Z masses of about 2 TeV and goes up to 20% for M Z 5 TeV. We observe that the impact of supersymmetric decays on the excluded M Z values runs from about 100 GeV (Z η ) to 200 GeV (Z ψ and Z N ), while the errors on the Z I dilepton rate in the UMSSM are too large to discriminate it from the non-supersymmetric case. Overall, Z bosons lighter than 4 TeV are still strongly disfavored by ATLAS data, regardless of the U (1) model.
In Fig. 2 (upper panel), we reexpress the same results by emphasizing the dependence of σB on the dilepton branching fraction, by superimposing the predictions of the different U (1) realizations, regardless of the actual θ E 6 mixing angle, and displaying the values of BR(Z → ll) by means of different colors. We find that the dilepton rate varies between 4% and 12%, and that the yielded exclusion masses are roughly between 4 and 4.5 TeV.
In the lower panel of Fig. 2, we present instead the distribution of the allowed Zboson masses as a function of the E 6 mixing angle, with the value of the g coupling for each scenario indicated by a color code. In order to determine the allowed regions, we first impose the experimental constraints in Table 4 and then the exclusion limits coming from the direct comparison with the ATLAS data in Fig. 1. The points ruled out by the ATLAS results are shown in grey. We observe, similarly to the findings of Ref. [27], that only |θ E 6 | values in the intervals [0, π/4] and [3/4π, π] can accommodate all the imposed experimental constraints. Outside of these regions, the U (1) charge of the extra singlet supermultiplet S is in fact close to zero so that either the SM-like Higgs boson or the Z boson, or even both, are predicted to be too light with respect to current data. In particular, Fig. 2 (lower panel) dictates that models U (1) χ and U (1) S are largely ruled out by the current data (see also the above discussion), while U (1) I is only marginally consistent. As a whole, after adding the recent ATLAS constraints [18] (the grey points), it turns out once again that scenarios exhibiting a Z boson lighter than 4 TeV can hardly ever be realized, the corresponding parameter-space regions getting more and more restricted.

Scenarios with Low-Scale Boundary Conditions
In this subsection, we focus on the second class of scenarios, wherein the input parameters, given in Eq. (2.17), are provided at the Z mass scale and where the U (1) coupling reads  Table 5, we learn that, for low-scale boundary conditions, g is substantially smaller. Therefore, the Z -production cross section is lower than for scenarios where boundary conditions are provided at the GUT scale M GUT .

2)
As a consequence, the inferred Z mass exclusion limits are reduced by about 200-300 GeV with respect to the high-scale unification case, as can be seen in Fig. 3, where the ATLAS limits are compared with the UMSSM predictions for U (1) ψ , U (1) η (upper panel), and U (1) N and U (1) I (lower panel) models. Since the g value is roughly the same as in the non-supersymmetric case, the overall impact of the inclusion of supersymmetric decays is similar to that found in the high-scale boundary framework, namely a reduction of the bounds on the Z boson mass by about 200 GeV. As observed for the other class of scenarios, the models with the highest impact of novel decay modes are the U (1) η and U (1) N ones, while the errors are too large to appreciate the effect of non-standard decays in σB for the U (1) I case. Our analysis then confirms the finding of Ref. [28], which compared UMSSM predictions in the low-scale unification framework with 8 TeV LHC limits and obtained an effect of similar magnitude on the excluded masses.
As for the high-scale unification case, we present in Fig. 4 (upper panel) the comparison of σB with the ATLAS data, scanning through the whole parameter space and displaying in different color codes the values of BR(Z → ll). Fig. 4 (lower panel) shows instead the correlations between the allowed M Z values and θ E 6 , accounting for both indirect constraints and direct ATLAS exclusion limits, the latter given by the grey-shaded area. The results in Fig. 4 are qualitatively similar to those presented in Fig. 2. However, as anticipated before, the g value is smaller, so that the ATLAS constraints on M Z are milder and values of M Z > ∼ 3.6 TeV are hence still allowed. Likewise, regarding specific U (1) models, U (1) χ and U (1) S are ruled out, while the other setups are still permitted and worth to be further explored.

Leptophobic Z Scenarios in UMSSM Models
The results presented in the previous section have shown that the inclusion of supersymmetric decays has a substantial effect on the Z searches and exclusion limits, but nevertheless the ATLAS bounds originating from the dilepton channel strongly constrain any phenomenologically viable UMSSM realization. Furthermore, the very fact that the Z boson has to be quite heavy impacts all sfermion masses through the U (1) D-terms, which may even lead to discarding some scenarios, such as U (1) χ , as yielding unphysical sfermion spectra. All LHC constraints studied so far can, however, be evaded by enforcing the Z boson to be leptophobic. In these scenarios, resonance searches in the dijet final state become the main probes of the new boson, Run II results for the top-antitop mode including the analysis of the full 2016 dataset being still not available. Dijet bounds are however much weaker, as described in Refs. [20,21].
Before discussing the phenomenology of leptophobic Z bosons within supersymmetry, in Fig. 5 we compare the CMS high-mass dijet yield from Ref. [21] with our predictions for σ(pp → Z ) × BR(Z → qq), obtained after scanning the UMSSM parameters as described in Table 3 and imposing the constraints of Table 4, for scenarios with high-scale (upper panel) and low-scale (lower panel) boundary conditions. As in the dilepton channel, the production cross section is calculated at NLO and the values of the dijet branching ratios are characterized by different color codes. For the sake of consistency with the experimental analysis, the σB rate is multiplied by an acceptance factor A 0.6 and the fraction of Z → tt events is not included in the calculation.
From Fig. 5, one learns that the computed σBA is always below the CMS exclusion limits in the range 2 TeV < M Z < 5 TeV at the 95% confidence level in both frameworks of coupling unification, once accounting for supersymmetric Z decays. One can, therefore, envisage than even much lighter Z bosons could be allowed by data when leptophobic UMSSM realizations, such as those introduced in Section 2.1, are considered.
Hereafter we focus on the second class of UMSSM scenarios, i.e., coupling unification at the M Z scale, and add to the list of free parameters in Eq. (2.17) the sine of the kinetic mixing angle sin χ, defined through Eq. (2.8), that we allow to vary in the [−1, 1] window. In principle, as thoroughly debated in Ref. [60], the kinetic mixing angle also affects the Dark Matter relic abundance, since the mass of the lightest supersymmetric particle (LSP), which in Ref. [60] can be either a right-handed neutrino or the lightest neutralino and in this paper isχ 0 1 , depends on the U (1) charges of the Higgs bosons, which have been Figure 5. Z production cross section multiplied by the dijet branching ratio and by the acceptance A 0.6, for the first (upper panel) and second (lower panel) class of scenarios investigated in this work. We compare NLO QCD theoretical predictions to the bounds obtained by the CMS collaboration [21] at the 1σ (green) and 2σ (yellow) level. The actual Z dijet branching ratio is indicated with the color code. modified according to Eq. (2.11) and are a function of sin χ 3 . Because of that, the authors of Ref. [60], besides applying the constraints due to collider physics, accounted for the upper bound on the relic density as well, relying on the Planck 2015 measurements [61]. The finding of Ref. [60] is that, although the mass of the LSP is indeed sensitive to sin χ and, e.g., a heavy Dark Matter candidate is favored by small | sin χ|, a value of the relic Figure 6. Correlations between the Z -boson branching ratio into a dilepton system and the θ E6 mixing angle featured by all points satisfying the constraints detailed in Section 2.2 and for UMSSM scenarios where the input parameters are fixed at the Z mass scale (second class of considered scenarios). The value of the sine of the kinetic mixing angle (sin χ) is indicated by the color code. density Ωh 2 ∼ 0.1, consistent with Ref. [61], can be achieved for any value of sin χ, and in particular for | sin χ| 0.3, corresponding to a leptophobic Z . In view of these results, we shall not impose further constraints, beyond those already discussed in the previous sections, and assume that any sin χ can possibly be consistent with the Dark Matter relic density, including the values which make the Z leptophobic.
In Fig. 6, we present the Z dilepton branching ratio, scanning the parameter space as presented in Section 2.2, in terms of the mixing angle θ E 6 and sin χ. In agreement with Eq. (2.11), we realize that values of sin χ around ±0.3 can lead to leptophobia whenever the E 6 mixing angle obeys the condition in Eq. (2.12) and the U (1) charges fulfill the relation Q E ≈Q L ≈ 0. In particular, the condition BR(Z → l + l − ) 0 can be achieved for −π < ∼ θ E 6 < ∼ − 3π/4, which includes the U (1) η model, and for π/8 < ∼ θ E 6 < ∼ π/4, hence in the neighborhood of U (1) ψ . The other U (1) symmetries are either ruled by the experimental data or, even in the most optimistic case, can hardly lead to dilepton rates below 5%.
Of course, these leptophobic scenarios cannot be constrained by standard Z -boson searches in dimuons or dielectrons at the LHC, and novel strategies must be designed. In the following, we propose a selection potentially allowing to observe leptophobic light Z bosons decaying through a supersymmetric cascade. As direct decays are forbidden, dilepton final states can arise from (Z -mediated) chargino-pair production and subsequent decays into a charged lepton and missing energy via an intermediate W boson, possibly off-shell, namelyχ ± 1 → (W ± → l ± ν l )χ 0 1 ,χ 0 1 being the lightest neutralino. However, for the points selected by our scan procedure, the off-shell contributions are typically either negligible (when the two-body decay channel is open) or not important enough to yield a sufficient number of signal events (when theχ ± 1 → W ±χ0 1 decay is closed). We, therefore, design an analysis assuming the presence of intermediate on-shell W bosons, targeting thus UMSSM scenarios where the mass difference between the lightest charginoχ ± 1 and the lightest neutralinoχ 0 1 is at least M W 80 GeV. The signal process consists of the resonant production of a chargino pair, followed by the decay of each chargino into a charged lepton and missing energy, We focus on two optimistic signal benchmarks that are currently not excluded by data and with different U (1) properties. Both scenarios exhibit a Z boson with a mass of about 2.5 TeV and charginos and neutralinos as light as possible, in order to maximize the branching ratios in Eq. (4.1), but with a mass splitting larger than M W , in such a way to allow the transitionχ ± 1 →χ 0 1 W ± with real W bosons. The first scenario, that we denote BM I, relies on a U (1) η symmetry, namely θ E 6 = −0.79π, since UMSSM scenarios based on this specific gauge symmetry can be made naturally leptophobic, as shown in Fig. 6. The second scenario, dubbed BM II, has instead a symmetry close to the U (1) ψ setup, but with a larger mixing angle, i.e. θ E 6 = 0.2π, so that a leptophobic Z boson can still be realized (see again Fig. 6).
The UMSSM parameters for the two points are quoted in Table 6, while Tables 7 and  8 contain the predicted masses for gluinos, squarks, sleptons, Higgses and gauginos in the reference points BM I and BM II, respectively. The branching ratios of the Z in such representative points are listed in Table 9, omitting rates which are below 1%. Table 6 shows that BM II features substantially larger values of tan β, µ eff and of the gaugino masses M 1 , M 2 , M 3 and M 4 , while M 0 and the trilinear couplings A 0 and A λ are smaller than in BM I. Comparing Tables 7 and 8, one learns that in BM I the squarks have masses between 3 and 4 TeV, while in BM II they are on average more than 1 TeV heavier. Charged sleptons in BM II are instead lighter than in BM I, unlike sneutrinos,   Table 7. Masses of gluino, squarks, sleptons, Higgs and gauginos for the UMSSM benchmark point BM I.q 1,2 ,l 1,2 andν 1,2 are mass eigenstates and differ from the gauge eigenstatesq L,R ,˜ L,R and ν L,R by virtue of the mass mixing contributions that are relevant especially in the stop case. All masses are in GeV.   Table 8. Same Table 7 but for the UMSSM benchmark point BM II.

Mg Md
whose masses vary between about 300 GeV and 1.7 TeV in BM I and between 660 GeV and 1.1 TeV in BM II. In the Higgs sector, with the exception of the SM-like h, all Higgs bosons have masses of a few TeV and are therefore too heavy to contribute to Z decays for both benchmarks.
In particular, as anticipated, the singlet-like neutral boson H has approximately the same mass as the Z , while H, A and H ± are roughly degenerate, with mass about 3.37 TeV in BM I and 5.24 TeV in BM II. As for gauginos, as anticipated, the two novel neutralinos χ 0 5 andχ 0 6 have masses similar to M Z , thus too high to be relevant for Z decays, while charginos and MSSM-like neutralinos are sufficiently light to possibly contribute to the Z width. Overall, the electroweakino spectrum is more compressed in the reference point BM I. The mass splitting betweenχ ± 1 andχ 0 1 is in fact slightly above M W in BM I, while it is substantially larger than M W , i.e. about 165 GeV, in the BM II framework. In both cases, the decayχ ± 1 → W ±χ0 1 can occur through on-shell W -bosons and has a branching fraction of almost 100%.
Concerning the Z branching ratios, Table 9 shows that the branching fraction of the Z boson decay into aχ + 1χ − 1 pair, entering in the process of Eq. (4.1), is of about 2% for the scenario BM I and 6% for the scenario BM II. BM I allows for substantial branching fractions into other combinations of chargino pairs, while both scenarios exhibit non-negligible rates into neutralino pairs, and the BM II scenario also includes decays into sneutrino pairs as well. The decay rates in pairs of the lightest neutralinos, possible candidates for Dark Matter, are instead suppressed in both reference points. As a whole, supersymmetric decays are responsible for 12% and 15% of the Z width in the representative points BM I and BM II, respectively.
Once our representative configurations are set, we carry out a full Monte Carlo event simulation at the LHC, for a center-of-mass energy √ s = 14 TeV. Hard-scattering signal events are generated with MadGraph5 aMC@NLO, the matrix elements being convoluted with the NLO set of NNPDF 2.3 parton densities. The production cross section is then σ(pp → Z ) 120 fb for both benchmarks. Parton showers and hadronization are simulated by means of the Pythia 8 program [62] (version 8.2.19), and the response of a typical LHC detector is modelled with the Delphes 3 package [63] (version 3.3.2), employing the Snowmass parameterization [64,65]. The resulting detector-level jets are reconstructed following the anti-k T algorithm [66] with a radius parameter R = 0.6, as implemented in the FastJet program (version 3.1.3) [67]. Moreover, we consider an average number of pile-up events of 140 and normalize our results to an integrated luminosity of 3000 fb −1 .
Regarding the backgrounds, we consider all processes leading to final states with two charged leptons and missing energy, such as vector-boson pairs V V , with V being a Wboson or a Z boson decaying leptonically. However, for the purpose of mimicking an actual experimental analysis, we account for processes yielding also jets which do not pass the acceptance cuts. Moreover, since our event simulation includes hadronization effects, we explore the possibility that background leptons originate from hadron decays as well. Overall, our backgrounds consist of single vector bosons (V ) or vector-boson pairs (V V ), possibly accompanied by jets, as well as tt and single-top events. In principle,  Table 9. Z decay rates for the benchmark points BM I (second column) and BM II (third column). Branching ratios below 1% are omitted.

-
even direct chargino production (pp → χ + 1 χ − 1 → l + l − + / E T ) should be considered as a background to the supersymmetric Z decays. Nevertheless, as pointed out in Ref. [26], the leptons produced in processes with direct charginos, unlike those coming from Z events, are typically pretty soft or collinear to the beams. It is therefore quite easy to suppress the pp →χ +  Moreover, in our selection strategy, we reject lepton candidates that are not at an invariant angular distance, in the transverse plane, of at least 0.4 from a jet, and only focus on muons that are cleaner objects than electrons, in particular for the pseudorapidity region considered in Eq. (4.2). We finally enforce the considered muons to be isolated, so that the activity in a cone of radius R = 0.4 centered on each muon contains at most 15% of the muon p T , We select events featuring two well-separated muons, since the two signal leptons l 1 and l 2 are expected to originate from two different supersymmetric cascade decays, by requiring N l = 2 and ∆R(l 1 , l 2 ) > 2.5 (4.5) and we veto the presence of jets, i.e.
Furthermore, the two signal leptons are expected to be produced from the decay of a heavy Z with a mass well above the TeV scale. We consequently impose the transverse momenta of the two leptons to fulfill p T (l 1 ) > 300 GeV and p T (l 2 ) > 200 GeV, (4.7) which are very efficient cuts to reduce the remaining SM background. We finally improve the sensitivity by requiring a large amount of missing energy, / E T > 100 GeV, (4.8) as could be expected for a signal topology where several neutrinos and neutralinos escape the detector invisibly.
The corresponding cutflows are shown in Table 10, which illustrates that, for the two benchmark scenarios under consideration, background rejection is sufficiently important for observing the signal despite the low selection efficiencies. For other possible benchmark choices (not considered in this work) featuring a heavier Z , the smaller production total rate is expected to be compensated by a larger efficiency of the two selection cuts restricting the transverse momenta of the two selected leptons.
Denoting the number of selected signal and background events by S and B ± σ B , we make use of two standard criteria, labelled as s and Z A , to define the LHC sensitivity to the leptophobic Z -boson signal, . (4.9) In Eq. (4.9), s is the significance as defined by the CMS Collaboration in Ref. [68] 4 , whereas the second method (Z A ) is known to be more suitable (and conservative) when the number of background events is small [70]. The conclusions are however very similar in both cases, 4 Following Ref. [69], the denominator of s sums in quadrature the intrinsic statistical fluctuation of the background √ B and the uncertainty in the background σB, thus obtaining s = S/ ( √ B) 2 + σ 2 B , leading to (4.9).
Step Requirements as can be seen from Table 10. For both significance definitions, we indeed find that the more compressed scenario BM I could lead to hints visible at the 3σ level, whilst the second scenario BM II is in principle observable at even more than 5σ. The largest LHC sensitivity to the latter scenario has a twofold origin. First, the Z -induced chargino-pair production cross section is larger by virtue of a greater BR(Z →χ + 1 χ − 1 ) branching ratio. Second, the heavier chargino mass typically induces harder leptons, the corresponding selection cuts being thus more efficient.
In the left panel of Fig. 7, we present the distribution in the transverse momentum of the leading muon l 1 after applying the first six cuts of Table 10. In the right panel of the figure, we in contrast show the transverse-momentum spectrum of the next-to-leading muon l 2 as resulting from the entire selection strategy. As for the p T (l 1 ) spectrum, all four considered backgrounds contribute at small p T , while above 100 GeV the only surviving SM events originate from the production of V V and tt pairs. The signal spectra are rather broad and lie below the backgrounds at low p T (l 1 ), whereas, for p T (l 1 ) > 300 GeV, both signals BM I and BM II start to be competitive with the background, yielding comparable numbers of events. For even larger transverse momenta, say p T (l 1 ) > 500 GeV, muons coming from supersymmetric decays of a leptophobic Z become dominant, especially in  Table 10 (left) and of the next-to-leading muon l 2 after applying all cuts (right) for both signal scenarios and the backgrounds. L dt = 3000 fb ∫ Figure 8. Left: missing transverse energy spectrum for the different components of the background and the two signal benchmarks. Right: cotransverse mass distributions for muon l 1 and invisible particles leading to missing energy (neutralinos and neutrinos). All histograms are obtained after applying all the acceptance cuts discussed in the paper. the reference point BM II. After all cuts are applied, the p T (l 2 ) distribution is explored (Fig. 7, right). All backgrounds are further suppressed and those due to single vectorboson and single-top production are negligible. The transverse momentum spectrum is thus substantial in the 200 GeV < p T (l 2 ) < 600 GeV range, with the BM II signal yielding the highest number of events through all p T range and BM I being also quite remarkable, especially for 200 GeV < p T < 400 GeV. Overall, Fig. 7 (right) shows that the cuts which we have applied are rather efficient to discriminate the leptons in leptophobic Z events from the Standard Model ones.
In Fig. 8 (left) we show the missing transverse energy, due to the lightest neutralinos χ 0 1 in our signal and to neutrinos in the backgrounds, after all cuts are imposed. The / E T spectra of our UMSSM benchmark scenarios are well above the backgrounds, once again limited to V V and tt pairs, through the whole / E T range. The BM II configuration, in particular, is capable of yielding a few events up to / E T 600 GeV, while, above 400 GeV, the backgrounds are basically all suppressed.
We have verified that any other transverse observable, such as the M T 2 or M CT variables defined in Refs. [71][72][73], are not useful for improving the considered selection strategy due to the too small mass difference between the lightest chargino and the lightest neutralino. The main features of the signal topology are in this case already captured by the requirements on the lepton transverse momenta and on the missing transverse energy. This is illustrated in Fig. 8 (right), where we present the cotransverse mass M CT distribution 5 for the leading muon l 1 and all particles contributing to the missing energy (lightest neutralinos and neutrinos). The M CT spectrum is qualitatively comparable to the / E T one. Both signals and backgrounds (V V and tt) peak at similar values, although the number of events generated by Z decays is always larger than for SM processes, and for M CT > 400 GeV only signal events survive. Designing an analysis with a possible extra cut on M CT would lead to a reduction in the significance, as both S and B would be affected in the same way. Such a new selection may, however, increase the sensitivity for spectra featuring larger mass gaps. In this work, we nevertheless choose to focus on the lighter UMSSM particle spectra that are still not excluded so far and thus more relevant for the near future.

Summary and Conclusions
Motivated by the latest ATLAS and CMS measurements which imposed improved lower bounds on the Z mass, we analyzed models with an additional U (1) gauge symmetry group arising from the breaking of E 6 supersymmetric GUT. We explored possible loopholes in the searches carried out at the LHC. In particular, we allowed the Z to decay into supersymmetric final states, such as gaugino pairs, and investigated scenarios where the Z is leptophobic. In fact, as the Z mass bounds are mostly determined by its decay into lepton pairs, the constraints would be relaxed in models in which direct leptonic decays are suppressed or even forbidden. We found that leptophobia can be achieved by accounting for the kinetic mixing between the two U (1) symmetries, parameterized by an angle χ, and that, among possible U (1) groups, the model U (1) η , while obeying all low energy conditions on the parameter space, is most favored to be leptophobic. Our analysis was then undertaken under two possible assumptions for scale unification, the gauge couplings being imposed to unify either at the GUT scale or at M Z . We investigated the mass bounds and decay patterns in both cases, as well as the prospects for seeing a Z signal above the background at the LHC, accounting for supersymmetry and leptophobia. 5 Given two particles of transverse energies ET,1 and ET,2 and transverse momenta pT,1 and pT,2, the cotransverse mass is defined as M 2 CT = (ET,1 + ET,2) 2 − ( pT,1 + pT,2) 2 [73].
Concerning supersymmetry, for both high-and low-scale unification, the rates of dilepton production are smaller once we include new decay modes, which translates into a reduction of the mass exclusion limits by about 200 GeV. As for dijets, we found an even larger impact of the inclusion of supersymmetric channels, so that the LHC constraints can be evaded. Within leptophobic scenarios, observing supersymmetric Z decays into charged leptons and missing energy would be most promising through a cascade from a primary decay into chargino pairs. We analyzed final-state signals from these intermediate states and suppressed the background by imposing a jet veto, in addition to requirements on the final-state leptons and missing energy. We chose two benchmark points in the parameter space, corresponding to different UMSSM realizations, and found that they both yield visible signals at the LHC, with a significance which varies from 3σ up to even 7σ, according to the criterion employed to estimate the LHC sensitivity.
Therefore, supersymmetric and possibly leptophobic Z decays are capable of giving detectable dilepton signals, which can be easily discriminated from the backgrounds and from non-supersymmetric Z events, so far employed to set the exclusion limits. Moreover, from the viewpoint of supersymmetry, Z bosons would be a promising source of new particles, such as the charginos and neutralinos investigated in this paper, which, unlike direct production in pp collisions, would feature additional kinematic constraints set by the high Z mass.
In summary, we believe that our study, accounting for Grand Unification Theories, supersymmetry and leptophobia altogether, should represent a useful guiding reference to explore a more general gauge structure than the Standard Model and address its incompleteness from perspectives that have not received so far proper consideration from the experimental collaborations. We demonstrated that investigating such scenarios is instead both worthwhile and feasible, as they are potentially capable of giving remarkable signals, especially in the high-luminosity phase of the LHC.