On New Physics Contributions to the Higgs Decay to $Z\gamma$

We explore models of new physics that can give rise to large (100% or more) enhancements to the rate of Higgs decay to $Z\gamma$ while still being consistent with other measurements. We show that this is impossible in simple models with one additional multiplet and also in well motivated models such as the MSSM and folded SUSY. We do find models with several multiplets that carry electroweak charge where such an enhancement is possible, but they require destructive interference effects. We also show that kinematic measurements in Higgs decay to four leptons can be sensitive to such models. Finally we explore the sensitivity of four lepton measurements to supersymmetric models and find that while the measurement is difficult with the high luminosity LHC, it may be possible with a future high energy hadron collider.

We explore models of new physics that can give rise to large (100% or more) enhancements to the rate of Higgs decay to Zγ while still being consistent with other measurements. We show that this is impossible in simple models with one additional multiplet and also in well motivated models such as the MSSM and folded SUSY. We do find models with several multiplets that carry electroweak charge where such an enhancement is possible, but they require destructive interference effects. We also show that kinematic measurements in Higgs decay to four leptons can be sensitive to such models. Finally we explore the sensitivity of four lepton measurements to supersymmetric models and find that while the measurement is difficult with the high luminosity LHC, it may be possible with a future high energy hadron collider.

I. INTRODUCTION
The Higgs boson's interactions with the electroweak gauge bosons, W , Z, γ are now well established. Decays to W W * [1][2][3][4], ZZ * [5,6] and γγ [7] have all been measured and are consistent with Standard Model (SM) expectations at the O(20%) level. Production of the Higgs in association with a W or Z [2,3] is also consistent with the SM. One decay mode that, as of yet, has not been detected is h → Zγ. The most recent searches for this channel from ATLAS [8] and CMS [9] place limits on the branching ratio relative to the SM prediction. The strongest limit comes from ATLAS which sets an upper upper bound on the h → Zγ signal strength at 3.6 times the Standard Model prediction.
In this work we explore the space of new physics (NP) models that can give large contributions to h → Zγ, close to the current experimental upper bound, without being in conflict with all other Higgs measurements. Due to precise measurements of the properties of the Higgs, and especially the Z and γ, we assume these fields are SM-like and that the NP contribution to h → Zγ arises at one loop due to new particles that couple to the Higgs and have electroweak quantum numbers. The SM contributions to h → γγ [34,35] and h → Zγ [10,11] are dominated by one-loop contributions (in contrast to h → ZZ, which is dominated by tree-level contributions), so oneloop NP can make O(1) modifications to the rate. Even in this case, however, the strong constraints on h → γγ make it difficult to engineer such models, and we show below that it is impossible to achieve this goal in models with a single new scalar multiplet. Therefore, if large deviations to h → Zγ are detected in the near future, this implies the existence of a complicated NP sector.
We will focus on NP scalars Φ with hypercharge and/or SU (2) L charge that contribute to h → Zγ. The case of other spins for the intermediate particles will not produce qualitatively different results. For a scalar Φ, regardless of its quantum numbers, one can always write renormalizable quartic interactions: where H is the SM Higgs doublet. 1 These operators induce one-loop contributions to h → γγ, h → Zγ, and h → ZZ of the type shown in Fig. 1. Trilinear operators written schematically as HΦH are also possible for certain quantum numbers of Φ, but these induce a vev for Φ which would give a significant contribution to the ρ parameter [36][37][38] that is excluded. With several multiplets with different quantum numbers, more complicated structures are possible: these will be discussed in Sec. III. New scalars with electroweak quantum numbers can be produced at hadron colliders such as the LHC and searched for directly, but the limits will depend strongly on their decay modes. The analysis here, however, is more model independent and only depends on the quantum numbers, the mass, and the coupling of the Higgs to the new states. At lepton colliders such as LEP, because of the relatively low background rates and fully known beam parameters, searches can be done in a more model independent manner [39][40][41][42][43]; these exclude masses below about 105 GeV. Given this lower limit, the Higgs cannot decay to an on-shell Φ, so we can expand the amplitudes from the diagrams in Fig. 1 in powers of m 2 h /4m 2 Φ . The leading term generated will be a CP -even dimension five operator: where F is the field strength tensor of the either the Z or γ.
Operators of the type given in Eq.
(2) also contribute to the Higgs decay to four leptons, the so-called golden channel, via diagrams of the form shown in Fig. 2. The leading contribution to this decay is the tree-level contribution mediated by ZZ * . Because this decay has four final state particles that can all be measured precisely, the rich final state kinematics can be used to gain significant information from each event. This means that one-loop contributions such as those from operators of the type in Eq. (2) can be probed [44][45][46][47] within the lifetime of the LHC. This has been used for various applications including the probing of exotic light states [48], measurement of the CP properties of the top Yukawa coupling [49], measurement of the ratio of the Higgs coupling to W W relative to ZZ [50], and the probing of various operators in SM effective field theory [51][52][53] as well as non-linear Higgs effects [54]. Even with the relatively small number of events already collected by the LHC, experiments can already use this data to place constraints on various scenarios [55].
Given models with large contributions to h → Zγ, we can estimate how well those models can be probed in h → 4 as a function of the number of events using the analysis techniques in [47,56]. We use a simple hypothesis testing procedure on a few benchmark points of these models and find that the high luminosity run of the LHC can reach approximately 2σ sensitivity to these models in h → 4 .
Rather than adding arbitrary new scalar multiplets, one can also ask if well motivated models can produce significant deviations in h → Zγ. Supersymmetric models are extremely well studied [57], and they contain scalar partners of the top quark, stops, which have large couplings to the Higgs and carry electroweak charges, so they can be scalars of the type shown in Fig. 1. The stops carry colour and, as such, contribute to production of the Higgs via gluon fusion in addition to the Higgs to diboson decays -measurements of these processes (especially h → γγ) can be used to place constraints on stop parameter space [58]. Here we update the constraints from [58] and show that these constraints imply that stops can only make small modifications to h → Zγ. We also find that h → 4 at the LHC will not be sensitive to these models, but there may be sensitivity at future higher energy hadron colliders.
Another well motivated model is folded SUSY [59], which also contains scalar top partners, F-stops, that have electroweak quantum numbers but do not carry colour. This makes direct bounds much weaker, and it also eliminates constraints from Higgs production via gluon fusion, making indirect bounds also weaker [58]. We also update the analysis on F-stops. While their contributions to h → Zγ can be larger than for ordinary stops, it still cannot be large, and the conclusions for the four lepton analysis is similar to that of stops.
The rest of the paper is organized as follows: Sec. II establishes the conventions used to describe new physics and the procedure used in analyzing the four-lepton processes, and Sec. III presents the construction and signals of models with large contributions to h → Zγ. Sec. IV presents an analysis of supersymmetric models and Sec. V wraps everything up.

II. HIGGS PHYSICS
Measurements of production and decay rates of the Higgs are all consistent with SM predictions. Therefore if there is new physics with electroweak charge that couples to the Higgs, it will be constrained by its contributions to the decay of Higgs to γγ at one loop. Since the leading SM contribution is also at one-loop, the constraints on such new physics at the weak scale will be strong. Similarly, new physics with colour charge that couples to the Higgs will contribute to the Higgs production via gluon fusion and can also place strong constraints.
Assuming the Higgs width is small compared to its mass, the cross section of a specific production mode i and decay to a specific final state V V can be parame- terized as [60]: Where Γ H is the total decay width of the Higgs, Γ V V is the partial decay width to V V , and σ i is the total cross section of (i → h → Anything) in the i production mode. Deviations from the SM can be parameterized through a set of coupling strength modifiers κ [7]. For a given production process or decay mode j, these are defined by: The SM has all κ j values equal to unity. The κ framework is not a full theory, and there are effects that it cannot account for that are detailed in [60]. The largest such effects are in processes where the Higgs is off-shell and not relevant for our analysis. We leave a more detailed accounting of the error budget due to this framework to future work. Generic NP contributions can be constrained by the limits placed by experimental measurements on Higgs boson coupling strength modifiers, κ j . In the models we are considering where the new physics dominantly contributes at one loop, the relevant coupling modifiers are those to photons, gluons, and h → Zγ: κ γ and κ g , and κ Zγ respectively. The kinematics of the measured processes do not change significantly in the presence of new physics above half the Higgs mass [7], so the κ framework is sufficient to describe deviations from the SM.
In Fig. 3 we show the constraints in the κ γ − κ g plane from [7] assuming all other couplings are SMlike. The best fit values resulting from their analysis are κ γ = 1.16 +0.14 −0.14 and κ g = 0.76 +0.17 −0.14 . On the other hand, the upper bounds on the Higgs decay to Zγ are, at 95% confidence level, 3.6 times the Standard Model prediction for the production cross-section times the branching ratio for pp → h → Zγ -this drops to 2.6 times the SM prediction if the production is set to the SM value. The best fit value for the signal yield is 2.0 +1.0 −0.9 [8], normalized to the SM. The future HL-LHC is expected to improve these constraints for the pp → h → Zγ process to 1.00 ± 0.23 the SM prediction [61]. In this work we will explore models that are not excluded by the κ γ − κ g analysis but that can have large contributions to h → Zγ, possibly saturating the current experimental constraint of 2.6 times the SM prediction.

A. h → 4 at One Loop
Despite the small rate, the h → 4 mode benefits from a rich kinematic structure and very high signal to background ratio. This, coupled with systematic uncertainties that are very different (and typically smaller) than with direct diboson measurements, make the four-lepton channel an important complementary way to examine the Higgs sector. Here, we follow the anaylsis of [47,56] in order to constrain these models using their contribution to h → 4 at one loop. The dominant contribution to h → 4 is via virtual gauge bosons with generic diagrams of the form of Fig. 2. The contributions to the hV V couplings can be described by an effective Lagrangian: Here, the leading order term that contributes to the hZZ coupling at tree level is given in L 0 . This term gifts the Z particle its mass and is generated via EWSB, where Z µ is the Z boson field. For the SM at tree level, we have A ZZ 1 = 2. Additionally, there are dimension-5 operators that can parameterize the one-loop contributions: where Z µν (F µν ) is the field strength of the Z (photon). These Lagrangian terms are only a subset of a more general description that would include other dimension five operators, but those operators are constrained to be too small to be relevant for this analysis [49,53]. Moving forward, we neglect the contributions of A ZZ 2 due to the lack of sensitivity [47,56,62] of future measurements to such modifications. It should also be noted that in the SM there are oneloop box and pentagon diagrams that contribute to the h → 4 process that are not captured by the parameterization of Eq. 7. These contributions are small, relative to the tree-level process, and are unaffected by the presence of any of the BSM models presented here. These calculations have been done [63,64], and we leave the integration of these results into our framework for future work. As such, we now turn our attention to the most relevant operators in the search for BSM: A γγ 2 and A Zγ 2 . Within the SM the A V γ 2 couplings will be generated primarily via a W boson loop along with (smaller) contributions from a top loop. Numerically, the SM values are A γγ 2 ≈ −0.008 and A Zγ 2 ≈ −0.014 [50]. 2 These form factors can be used to compute on-shell Higgs decay to γγ or Zγ. These SM one-loop contributions have been explored in the literature for both h → Zγ [10,11] and h → γγ [34,35]. We have assumed NP contributions are CP -even because of strong constraints from flavour and CP violating observables [66]. 3 Given this parameterization, we express the matrix element for the h → 4l process in the form, where V = Z, γ, and P µν (V i ) are the propagators of the vector bosons. It should be noted that the first line is where any potential NP that we could hope to uncover would reside. The second line is simply vector bosons propagating and decaying to leptons -processes that are well understood and well measured. The h → V 1 V 2 2 Note that [65] gives a different sign for A Zγ 2 . 3 For analyses of CP violation in h → Zγ see [67][68][69].
matrix element can be parameterized as follows: with i = ZZ, Zγ, γγ and k 1 , k 2 representing the four momenta of the intermediate vector bosons (or, equivalently, lepton pairs). The form factors, C i n , are Lorentz invariant and encode the momentum dependence. They have the form, with f i (m 2 h /m 2 X , k 2 1 , k 2 2 ) representing the loop function for NP particle X coupled to the Higgs with coupling h X . Previous studies have indicated that dependence on the invariant mass parameters k 2 i is quite weak [49,70]. In order to examine this claim, the form factors C i n were expressed as a Taylor Series and expanded around the pole masses of the intermediate vector bosons. The parameter space examined here can indeed be probed taking k 2 i = m 2 for the relevant boson, i.e. treating the vector as on-shell for the form factor. This limit is used for determining the exclusion bounds from ATLAS measurements of Higgs to diphotons.
All the NP scenarios we will explore can then be encoded in the value of the effective field theory (EFT) operators A γγ 2 and A Zγ 2 for all models (plus the gluon form factor A gg 2 for coloured new scalars such as stops). They are computed via the diagrams in Fig. 1  In the Higgs rest frame, the useful kinematic variables for this decay are: • Φ: The decay angle between the decay planes of the intermediate bosons in the rest frame of the Higgs.
• θ 1 : The angle between the lepton coming from the decay of V 1 and the momentum of V 2 in the V 1 rest frame.
• M i : The invariant mass of the lepton pair produced by the decay of the i th boson.  [71]. Since gluon fusion is the primary production mechanism for the Higgs at the LHC [72] and the variables examined during our analysis are in the rest frame of the Higgs and thus only sensitive to the decay process, including other production modes would not change our analysis. Additionally, the very high signal to background ratio of the four lepton channel and the easy identification of leptons within LHC detectors makes both a full simulation of backgrounds and an examination of detector effects subleading for the purposes of this work [62]. Measurements (see for example [73]) in this channel confirm that the signal to background ratio is very large in the Higgs invariant mass window.
We now determine the number of events needed to distinguish two different hypotheses (SM vs. NP) at a given confidence level. This is done with a likelihood analysis of the MC generated NP events compared to MC events assuming purely SM physics. A standard unbinned likelihood analysis is used (see [45,74] for further details) where the computed differential cross-section can be normalized and used as a probability over the 4l kinematic variables: where Y is the set of kinematic variables of an event, A is the given set of effective coupling strengths for the model being examined, and M is the matrix element from Eq. 8. The integral is over the fiducial acceptance region for 4 events. If we have N events, this allows us to compute a likelihood: If we have two different scenarios (typically a SM scenario and a NP scenario), we can calculate a likelihood ratio and construct a hypothesis test statistic, Λ = 2 log (L(A 1 )/L(A 2 )) .
This test statistic is calculated on MC data generated assuming each of the two underlying hypotheses. This creates two distributions of the test statistic Λ, one assuming scenario 1 is true, the other, scenario 2. The separation between these distributions is a measure of how easily distinguishable the two hypotheses are. Explicitly, if we call the distribution with the smaller average test statistic f and the other distribution g then there exists a value Λ 0 such that, Because Λ 0 is the point at which the distributions are indistinguishable, the value on either side of Eq. 13 is the one-sided Gaussian probability of the observed events of one scenario excluding the alternative. This can be converted into a more standard σ value that is a function of the number of signal events N . In other words, we can get the expected statistical significance as a function of the number of events for any two scenarios.

III. MODELS WITH LARGE h → Zγ CONTRIBUTIONS
Here, we explore models with new electroweak multiplets that can give large contributions to h → Zγ while not being in conflict with h → γγ.

A. New Scalar Multiplets
We begin with a simple scenario: a single scalar multiplet that is charged under SU(2) × U(1), but is a singlet under SU(3) c . At the moment, we remain agnostic to the size of the multiplet; we do note, however, that the size of the multiplet is bounded at 8 (isospin of 7/2) due to perturbative unitary constraints [75]. Finally, in order to stay away from any LEP II bounds [39][40][41][42][43], we require our particles have a mass 100 GeV.
The NP contributions to our effective operators A Zγ 2 and A γγ 2 have the form where h X is the coupling of the new particles to the Higgs, g X and z X are the coupling of the new particles to photons and Z bosons, respectively, and l γγ and l Zγ are the loop functions that depend on the mass of the NP particles given in App. A. The sum is over all states in the multiplet. It should be noted that, unless there is mass mixing, h X and the particle mass is the same for all members of the multiplet. The couplings g X and z X , on the other hand, depend on the T 3 value for the particle and are thus different for every multiplet member. Additionally, for NP particles with mass 100 GeV, the loop functions become approximately equal, l γγ [m] ≈ l Zγ [m], so the relative contribution sizes of the new particle multiplet to the diphoton vs. Zγ is controlled by the relative size of with the sum being over all particles in the multiplet. The relevant terms in the Lagrangian for a single scalar multiplet X are given by with D µ being the typical covariant derivative that gives rise to the coupling of X to gauge bosons. The coupling b gives rise to the Higgs coupling to X pairs after electroweak symmetry breaking, and we ignore X selfinteractions as they do not affect the analysis. We assume m 2 X > 0, as bounds on new electroweak scalars with vacuum expectation values are very strong [36][37][38] (unless X has the same quantum numbers as the Higgs). For the same reason, we do not include HXH terms which are allowed for certain quantum numbers of X. Finally, we remain agnostic as to which decay modes are present for the new particles; this allows us to ignore direct search limits and focus wholly on Higgs decays.
After a rotation into the mass basis for the gauge bosons, the aforementioned couplings for the new scalars are: where g and g are the gauge couplings for SU (2) L × U (1) Y and s W and c W are the sine and cosine of the weak mixing angle. The typical convention relating electric charge of the scalars to the isospin and hypercharge holds: Requiring integer charge and half-integer (integer) values for weak-isospin restrict hypercharge to half-integer (integer) values.

B. The Failure Of One Additional Multiplet
We now demonstrate that the inclusion of a single multiplet cannot give rise to large contributions to h → Zγ. The relevant products of couplings to gauge bosons are given by Summing over T 3 gives effective couplings (defined in Eq. (15)) where k is the size of the multiplet and Since g 2 cos 2 θ W > g 2 cos θ W sin θ W , g ef f γγ will always increase in strength faster than g ef f Zγ as the absolute value of the multiplet's hypercharge becomes larger. So, to maximize the relative size of the Zγ contributions, the magnitude of the multiplet's hypercharge must be kept as small as possible Therefore, the best case scenario is an odd multiplet with Y = 0 that gives The even multiplets give weaker results, though as k becomes larger, this ratio approaches the same value as the odd multiplets. When we consider the fact that in the SM the Zγ effective operator is larger than the γγ one, this shows the modification of to the h → Zγ and h → γγ rates are at best comparable. As such, the large h → Zγ parameter space cannot be probed by such simple models and we must add another degree of complexity through the addition of more multiplets. These scenarios can result in large h → Zγ rates through the cancellation of h → γγ contributions and/or with mass mixing effects such that the mass eigenstates with comparatively small γ couplings have low masses and those with larger electric charges have larger masses. The mixing scenario turns out to be less compelling: even if one were to obtain mass mixing such that the member of the multiplet with the largest T 3 value (e.g. +1/2 for a doublet) that has the largest contribution to h → γγ relative to h → Zγ is effectively decoupled, the improvements to Eq. (21) are minor. Even with more complicated scenarios, it is extremely difficult to construct a model with a viable scalar potential that achieves the required mass splitting. With this in mind, for the rest of this section we explore cancellation scenarios.

C. The Singlet-Triplet Model
The road map of where to go next comes from analyzing this question from the point of view of the unbroken electroweak effective field theory. In the EFT prior to spontaneous symmetry breaking [76][77][78] (typically called SMEFT) there are three dimension-6 terms relevant to Higgs to electroweak diboson decays [79]: where W µνa and B µν are the field strength tensors of SU (2) L and U (1) Y respectively, and τ a are Pauli matrices. The introduction of scalar multiplets can generate contributions to the Wilson coefficients Q of either sign: this can lead to multiplets that either enhance, dampen, or even nullify each other's effects.
We can get expressions for our dimension-5 EFT in the Higgs basis defined in Eq. (7) in terms of our unbroken operators: In the 3-dimensional space of the of Q W , Q B , and Q W B , there is a plane of values that lead to no contribution to A γγ 2 -yet the vast majority of this plane (everything outside one line) does offer up contributions to A Zγ 2 . The task, then, becomes finding an appropriate set of multiplets that can simultaneously enhance the rate of h → Zγ and leave h → γγ close to its SM value.
A model that features these properties is one where we introduce two NP scalar multiplets: a singlet state, S, with Y = 1 and a triplet state, T , with Y = 0. Both multiplets are singlets under SU(3) c . At leading order, integrating out S generates Q B and integrating out T gives Q W . Therefore, by choosing couplings appropriately, we can be on the line in the Q B − Q W plane where the contributions to h → γγ (approximately) vanish. As we will consider relatively light states, we will use the EFT as only a guide and do our analysis in the full theory.
The relevant NP Lagrangian terms are: The contributions to the Higgs to diboson decays are simply those from Eq. (14) generalized to allow for multiple particles, where h S and h T are equal to b v and d v, respectively; while g ef f T γγ and g ef f T Zγ are calculated through Eq. (19). Just as in the single NP multiplet case, the particle masses we consider are greater than 100 GeV and, as a result, the loop functions for γγ and Zγ are nearly identical.
Cancellation of the NP contributions to A γγ 2 requires either h S or h T (or, in terms of the Lagrangian, b or d) to be negative. For concreteness we select h S to be negative. The scalar potential must be bounded from below in order to maintain its stability: this can be checked by ensuring that the quartic terms of the scalar potential have a limit V s ≥ 0 in all directions (any direction where V s = 0 also requires the limit of the quadratic terms in that direction to be ≥ 0). This requires the quartic couplings a and c to be positive, and we have also taken d > 0. As it does not affect our analysis, we can also take e > 0.
The only remaining constraint is b < 0 potentially leading to an unbounded direction. This analysis is nearly identical to those in 2 Higgs Doublet Models [80]. Setting H † H = r cos θ and S † S = r sin θ, followed by taking the large r limit gives r 2 a sin 2 θ + λ cos 2 θ + b cos θ sin θ = r 2 f (θ) (27) where λ is the SM quartic Higgs coupling. Ensuring Eq. (27) is positive can be done by requiring f (θ) ≥ 0 at its smallest point. Applying this requirement enforces a minimum size on the negative coupling b: The singlet-triplet model is representative of a much larger class of possible models -not a unique solution to a large Zγ signal. This construction can apply to two multiplets of any sizes, although this particular choice of charges is simpler than the generic one because there are no allowed mass mixing terms. Within this model, we will work with the following benchmark parameter point: This point evades constraints from LEP [39][40][41][42][43], has reasonable coupling sizes, scalar potential stability, and is within the 2σ limits on h → γγ [7]. Full cancellation of the NP diphoton decay contributions is possiblehowever to avoid tuning parameters we are content with values that respect current bounds. This point also saturates ATLAS [8] upper limit on h → Zγ. We take this as a representative of the relatively complicated models required to generate a large anomaly in h → Zγ while being consistent with other constraints.

D. Four-Lepton Sensitivity
Here we calculate the number of h → 4 events required to probe the above benchmark point in the singlettriplet model using the analysis described in section II B. From the loop diagrams in Fig. 1, we can determine the contributions to the effective operators A γγ 2 and A Zγ 2 : The modification of the branching ratio (or event rate) is given (1 + A(NP)/A(SM)) 2 . In terms of κ γ and κ Zγ : Applying the 4 analysis, we obtain the results shown in Fig. 4. Specifically, we generated 1 million MC sample pp → h → 4l events for both the SM scenario and our singlet-triplet benchmark, and the liklihood was constructed using multiple pseudo-experiments with varying number of events N . In Fig. 4, the points in the left figure show the discrimination power between to the SM and NP models presented in terms of Gaussian σ as a function of N , the number of events. We also impose a fit curve of the form x 0 + x 1 √ N as √ N growth is expected for large N with large signal to background, as is the case in h → 4 .
The curve from the left figure is extrapolated in the right part of Fig. 4 in order to estimate the required number of events to distinguish these two scenarios. With ∼ 8000 events, we can rule out the parameter point at 2σ and, more generally, we can begin to probe the large Zγ model parameter space. Using 139 fb −1 of data, ATLAS has measured the fiducial cross section of the H → ZZ * → 4l to be σ f id = 3.28 ± 0.32 fb in good agreement with the SM prediction σ f id = 3.41 ± 0.18 fb [81]. With the high luminosity (HL) LHC expected to collect upwards of 3000 fb −1 over its lifetime [82], the LHC should be able to begin seriously probing the large Zγ parameter space by the end of its run.
At a future high energy hadron collider the prospects improve further. At 100 TeV center of mass, the expected gluon fusion production cross section is σ = 808.23 +44.53 −56.95 pb [83] with additional contributions of ∼ 10% total from vector boson fusion and associated production [84]. Taking the h → 4l branching ratio to be 2.796×10 −4 [85] and assuming an integrated luminosity of 10 ab −1 gives ∼ 2.5 million signal events -enough to potentially probe the relevant parameter space of this model. With this level of precision, other effects not considered here such as backgrounds, detector effects, and higher order corrections may become important.

IV. SUPERSYMMETRIC MODELS
We now turn to study supersymmetric models [57] to see if large effects in h → Zγ can appear in well motivated models. Using the analysis techniques of Sec. III, we examine scalar top partners which, in the MSSM, generically have the largest coupling to the Higgs of new supersymmetric states. Unsurprisingly, the additional restrictions coming from the extra structure present in SUSY models leads to much smaller allowed contributions to h → Zγ than with the general multiplet scenarios. In addition, for SUSY models we find that the deviations in the branching ratio of h → Zγ are always smaller than the deviations in h → γγ.
In addition to considering stops in the MSSM, we also consider F-stops, scalar top partners in folded SUSY models [59]. This class of models stabilizes the weak scale against radiative corrections to around 5 TeV. The scalar partners in these models are not charged under SU(3) c , but have the same electroweak quantum numbers as stops allowing them to couple to the Higgs in the same way. This results in F-SUSY contributions to the Higgs decays to two photons and Zγ, but no modification to the gg → h amplitude (although there is a new decay to hidden gluons, h → g h g h ). There are many versions of F-SUSY to chose from; here we take the simplest approach and imagine a scenario identical to the MSSM stop sector sans stop-gluon couplings [58]. We examine the current restrictions on the stop parameter space for both models and then complete a four lepton analysis to estimate the required number of Higgs events to probe these scenarios.
The (F)-stop parameter space is dominantly controlled by three parameters: the two mass eigenvalues m 1 , m 2 , and the mixing angle between the gauge and mass eigenstates, θ t (we use the conventions in [86]). There is some weak dependence on tan β -the couplings change less than 10% over any possible β value and, if we restrict this parameter to its normally accepted values, 3 ≤ tan β ≤ 50, the modification becomes O(1%). As such, we take tan β = 10 for simplicity.
Using standard formulas in the literature we can compute the contributions of stops to gg → h and h → γγ and apply the ATLAS constraints shown in Fig. 3. We ignore direct bounds on (F)-stops in order to keep our constraints model independent. As this is a three dimensional parameter space, we present our results in Fig. 5 in two dimensional slices of the two mass eigenstates with the mixing angle fixed at θ t = 0 (left) and the lighter mass eigenstate m 1 vs. the mixing angle θ t with m 2 = 1 TeV. We note that these results are leading order, and higher order corrections place an uncertainty on these bounds of O(10%). This is an update of the analysis in [58] using significantly more data. Scanning the parameter space indicates that MSSM stops below ∼ 140 GeV are essentially excluded. Additionally, both stops having a mass under ∼ 200 GeV is excluded. In general, the weakest parameter space constraints occur when either mixing is maximized or the stops have a very large mass splitting and are effectively decoupled.
We can now compute the contribution to the CP -even form factor A 2 from Eq. (9). 4 By scanning over the m 1 , m 2 , and θ t parameter space and subjecting the results to the previously mentioned ATLAS constraints, the maximal modification of said form factor (and thus matrix element) was found to be: This is beyond the precision expected on the coupling κ Zγ from the HL-LHC of O(10%) [61], but deviations from the SM in this channel may be visible at a 100 TeV collider with an expected precision of O(1%) [88], but this would require reducing the uncertainty in the SM prediction which is currently O(5%) level [85]. If this is the theory of nature, we expect to see deviation from the SM in h → γγ 5 long before we see one in h → Zγ. For F-stops, the procedure is the same, but we can ignore constraints from gluon fusion production of the Higgs. This makes measurements less constricting leading to the results presented with the dashed lines in Fig. 5. The notably weaker constraints lead to stops only below ∼ 100 GeV being universally excluded and requiring at least one stop above ∼ 150 GeV. Other than the strength of the bounds, the overall takeaways remain the same as with MSSM stops: the weakest constraints on the lighter stop occur when either mixing angle has values of θ t = 0, ± π 2 or the stops have a very large mass 4 CP violation can exist in the MSSM but the constraints are strong [87]. 5 Note that these deviations are significantly larger than the uncertainty on the SM prediction which is of O(1%) [85].
splitting. The A 2 values in the F-stop allows greater deviations: but we see that the deviations in h → Zγ are much smaller than the singlet triplet model shown in Eq. (30). This is just on the edge of the sensitivity of the HL-LHC in on-shell h → Zγ [61], but deviations would first be seen in h → γγ.

A. Four-Lepton Analysis
The four-lepton analysis of the SUSY models follows the one in Sec. III. For stops, we use the following benchmark point: where we have introduced the gluon form factor A gg 2 analogous to those of the electroweak gauge bosons. Simi- As in the large Zγ case, these coupling modifications were used in generating MC events that were then processed into likelihoods, test statistic distributions, and, finally, statistical confidence as a function of number of signal events σ(N ). The results for both stop and Fstop scenarios are shown in Fig. 6. Unsurprisingly, the much smaller impact of the intermediate scalar particles within the loop decays makes both the stop and F-stop much more difficult to probe than the large Zγ models. With the HL-LHC's integrated luminosity topping out around 3000 fb −1 after its final run [82], the required ∼ 150 thousand events to be sensitive F-stops or ∼ 270 thousand events required to explore stops at 2σ are not reachable with the LHC. These numbers are, however, potentially in reach of a future 100 TeV collider -having an integrated luminosity on the order 10 ab −1 gives ∼ 2.5 million signal events which can potentially probe the SUSY parameter space.

V. CONCLUSION
We have explored the possibility of new physics models with significant enhancements to the h → Zγ decay. We have demonstrated that simple models cannot give rise to large enhancements while still being consistent with other data, particularly measurements of h → γγ. Models with multiple multiplets featuring auspicious cancellations could produce such a signal. For models that do feature large contributions to h → Zγ, such as a model with a singlet and a triplet explored in section III C, kinematic analysis of h → 4 will be able to probe these models with the data from the high-luminosity LHC.
We also explored more motivated models that can solve the hierarchy problem, focusing in particular on stop contributions in supersymmetry and colourless top partners (F-stops) in folded SUSY. We used current Higgs data to constrain the parameter space of those states updating the analysis of [58], and determined these models cannot give measurable deviations to h → Zγ. Kinematic analysis of h → 4 can only have sensitivity to these models with significantly more data than the LHC will have although this may be possible with a future 100 TeV hadron collider.
Although this work focused exclusively on the new physics contributions of scalars, an analagous analysis with fermionic multiplets would not change the qualitative conclusions: the loop factors and Higgs couplings are slightly different, but the photon and Z couplings still arise from the covariant derivative and models with large contributions to h → Zγ can only be produced through fortuitous cancellation by multiple new physics multiplets. CP violating couplings were also not examined in this work, but the constraints on Higgs couplings due to EDM measurements [66] indicate that such couplings would have to be tiny -much too small to achieve the large Zγ contributions that we examined.
What this ultimately means is that, in general, the h → Zγ channel is unlikely to be a place where new physics is discovered: simple or motivated models tend to lead to contributions smaller than the much more sensitive diphoton channel. As such, a discovery of a significant NP contribution to this channel is a strong indication that interference effects are present in the NP sector and points towards non-minimal models like those presented in this paper. Finally, the discriminating power of a 100 TeV collider, particularly in h → 4 , is undeniable: where the LHC cannot even probe the most favourable parameter points of the MSSM, a future collider will collect significantly more events, potentially giving us a much clearer handle on the nature of our microscopic world. For the Zγ expression, the case presented here de-scribes loops with only one type of new particle; loops where the legs are different particles lead to longer expressions we do not present here. In the case of the singlet-triplet model presented in Sec. III C, this expres-sion holds -however for the SUSY models in Sec. IV the more general expression is required due to the possibility of the presence of both stops.