Precision from the diphoton Zh channel at FCC-hh

The future 100 TeV FCC-hh hadron collider will give access to rare but clean final states which are out of reach of the HL-LHC. One such process is the $Zh$ production channel in the $(\nu\bar{\nu} / \ell^{+}\ell^{-})\gamma\gamma$ final states. We study the sensitivity of this channel to the $\mathcal{O}_{\varphi q}^{(1)}$, $\mathcal{O}_{\varphi q}^{(3)}$, $\mathcal{O}_{\varphi u}$, and $\mathcal{O}_{\varphi d}$ SMEFT operators, which parametrize deviations of the $W$ and $Z$ couplings to quarks, or, equivalently, anomalous trilinear gauge couplings (aTGC). While our analysis shows that good sensitivity is only achievable for $\mathcal{O}_{\varphi q}^{(3)}$, we demonstrate that binning in the $Zh$ rapidity has the potential to improve the reach on $\mathcal{O}_{\varphi q}^{(1)}$. Our estimated bounds are one order of magnitude better than projections at HL-LHC and is better than global fits at future lepton colliders. The sensitivity to $\mathcal{O}_{\varphi q}^{(3)}$ is competitive with other channels that could probe the same operator at FCC-hh. Therefore, combining the different diboson channels sizeably improves the bound on $\mathcal{O}_{\varphi q}^{(3)}$, reaching a precision of $|\delta g_{1z}| \lesssim 2 \times 10^{-4}$ on the deviations in the $ZWW$ interactions.

Abstract: The future 100 TeV FCC-hh hadron collider will give access to rare but clean final states which are out of reach of the HL-LHC. One such process is the Zh production channel in the (νν/ + − )γγ final states. We study the sensitivity of this channel to the O (1) ϕq , O (3) ϕq , O ϕu , and O ϕd SMEFT operators, which parametrize deviations of the W and Z couplings to quarks, or, equivalently, anomalous trilinear gauge couplings (aTGC). While our analysis shows that good sensitivity is only achievable for O (3) ϕq , we demonstrate that binning in the Zh rapidity has the potential to improve the reach on O (1) ϕq . Our estimated bounds are one order of magnitude better than projections at HL-LHC and is better than global fits at future lepton colliders. The sensitivity to O (3) ϕq is competitive with other channels that could probe the same operator at FCC-hh. Therefore, combining the different diboson channels sizeably improves the bound on O

Introduction
In the coming decades, future high-energy colliders will provide the next step in the study of the Standard Model (SM) and in the exploration of new physics. An accurate preliminary assessment of their physics potential is, therefore, necessary to select the most promising machines and set up a tailored and effective experimental strategy.
Among the various proposals for future accelerators, high-energy hadron colliders play a dominant role due to the broad range of new-physics probes they allow for. Although their primary target is indubitably the direct search for new massive particles, they also offer a complementary way to look for new physics through indirect probes. The latter strategy already proved successful at the LHC, where indirect electroweak (EW) probes were able in some cases to surpass the accuracy achieved at LEP [1,2].
The success of the indirect search strategy is based on the fact that new-physics corrections, below the threshold for direct production, tend to grow with the energy of the process. Hadron colliders have the advantage of probing a large range of energies, thus allowing us to test the high-energy tails of kinematic distributions, where new physics is more easily accessible. The exploitation of selected "clean" processes with small statistical and systematic uncertainty is key for this kind of precision studies, especially when they target the EW dynamics [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16].
One of the main limiting factors in carrying on the EW precision program at the LHC is the limited amount of events in many clean channels, which forces one to focus only on final states with high cross-section. Some of the best probes for new physics come from processes involving gauge bosons and/or the Higgs, most notably the di-boson production channels [2][3][4][5][6][7][8][9][10][11][15][16][17][18][19][20][21][22][23]. These processes are sensitive to many operators that describe the couplings of the SM gauge bosons and the high-energy dynamics of the Higgs. At the LHC only channels with large decay branching ratios are accessible. For the Higgs, in particular, this means that the analysis is basically limited to the h → bb decay mode, with the necessity to cope with large backgrounds and a relatively complicated final state. For instance reconstructing an energetic Higgs in the bb channel requires the use of sophisticated boosted-Higgs analysis techniques.
Future high-energy hadron colliders offer a big advantage from this point of view, since they provide significantly larger cross-sections and much higher integrated luminosity. This allows one to broaden the set of decay channels suitable for high-energy precision probes. An intriguing possibility is to consider the associated production of a Higgs with a gauge boson, in which the Higgs decays in an easily reconstructable final state. For instance the h → γγ decay mode, especially if combined with a leptonic decay of the EW boson, can offer a very clean signature with virtually no background and small systematic uncertainties.
In a previous paper by a subset of the authors, Ref. [15], we already started an analysis of this kind of processes focusing on W h production at a future 100 TeV proton-proton collider (FCC-hh). We considered the final state in which the Higgs decays into a pair of photons, while the W boson decays to light leptons (electrons and muons). We found that good sensitivity to new-physics effects that modify the coupling of the W -boson to quarks can be achieved. Interestingly, the expected accuracy can surpass the one achievable at the end of the LHC program by one order of magnitude and is competitive with EW precision measurements at future high-energy lepton colliders [15] (for instance FCC-ee).
In this paper, we continue the study of the associated Higgs production channels, turning our attention to the Zh production process. As in Ref. [15], we focus on the channel in which the Higgs decays into a pair of photons. For the Z boson, we consider two possible decay channels, namely the one into a pair of light charged leptons and the one into neutrinos. As we will show, both channels offer very distinctive signatures that can be used to make the analysis almost background free.
To parametrize the new-physics effects, we adopt the Effective Field Theory (EFT) approach, assuming that the beyond-the-SM (BSM) dynamics has an intrinsic scale sufficiently higher than the range of energy spanned by the Zh events. As we will show, in the case of FCC-hh, this amounts to the requirement that the EFT cut-off is 2 TeV (see With respect to W h production, the Zh process is sensitive to a wider set of energygrowing new-physics effects. In fact it can be used to test the same operator, O (3) ϕq , that gives the leading corrections in W h. Moreover it gives access to three additional operators, O (1) ϕq , O ϕu and O ϕd , that encode deviations in the Z-boson coupling to left-handed and right-handed quarks. All these operators give rise to corrections that grow quadratically with the center-of-mass energy of the process, and are thus more visible in the high-energy tail of the kinematic distribution. This paper is structured as follows. In Section 2, we discuss the properties of the Zh production channel and the leading new-physics effects it is sensitive to. In Section 3, we discuss the details of our analysis. In particular, we present the features of the signal and background processes and the cuts we exploited to single out the new-physics effects from the backgrounds. In Section 4, the results of our analysis are collected, whereas the conclusions of our work and possible future research directions are discussed in Section 5. We collect in Appendices A and B some explicit formulae and additional details about the analysis that were not included in the main text.

Leading new-physics contributions
We parametrize new-physics effects via the SMEFT formalism, focusing on the leading BSM effects due to dimension-6 operators. Only four operators lead to contributions to the pp → Zh amplitudes that grow quadratically with the partonic center-of-mass energy [2]. In the Warsaw basis [24], they are given by and σ a are the Pauli matrices. We define the Wilson coefficients of the above operators to be dimensionless so that the effective Lagrangian can be written as where Λ is the EFT scale and the index a runs over the operators in Eqs. (2.1)-(2.4).
To report all the numerical results we will fix the EFT scale to the conventional value Λ = 1 TeV.
While, in principle, each operator also carries a generation index, we focus on the flavoruniversal scenario in which the operators couple diagonally and with the same strength to all quarks. This assumption is physically motivated by the reduction of flavor constraints, in particular the stringent ones coming from flavor-changing processes. As we will argue later on, the bounds from our analysis mostly come from the couplings to first-generation quarks. Therefore, our results remain valid to good approximation if the operators couple only (or mainly) to the up and down quarks.
We also neglect any modification to the Higgs and Z-boson branching ratios, specifically h → γγ, Z → + − , Z → τ + τ − , and Z → νν. This approximation is justified by the fact that, by the end of the FCC-hh program, these branching ratios are expected to be known with a precision comparable to (or below) the luminosity uncertainty. The branching ratios of the Z boson are currently measured with per mille precision [25]. 1 On the other hand, the bound on the deviations in the effective hγγ coupling is expected to be at most 2% by the end of HL-LHC and below 0.3% by the end of FCC-ee+FCC-hh [26]. We also note that, due to the energy enhancement, the sensitivity to the effective operators (2.1)-(2.4) mainly comes from the differential distribution in the process energy. The overall cross-section, which is dominated by the low-energy phase-space region, is only weakly affected by new-physics contributions and plays a marginal role in the fit.
We are mainly interested in studying small deviations from the SM, i.e., the case where the SM contribution to the amplitude dominates over the BSM one. In such a case, the main deviations from the SM predictions are expected to come from the interference between the SM and the new-physics contributions, since the square of the BSM amplitude is of higher order in the EFT expansion and therefore more suppressed. Optimizing the sensitivity to the interference terms is therefore a crucial step in the analysis. This generic picture might fail in specific situations, however. As we will see in the following, in some cases the interference terms might be small or altogether absent due to structural or accidental cancellations. When this happens, new-physics effects come dominantly from squared BSM amplitudes which leads to two additional complications. First, new-physics effects become measurable only for relatively large values of the Wilson coefficients. And, second, additional contributions arising from dimension-8 operators, which we do not consider in our analysis, could play an important role and make the bounds more model dependent [27].

Structure of the interference terms
We start by discussing the high-energy behavior of the leading-order (LO) pp → Zh helicity amplitudes. The scaling with the center-of-mass energy of the process,ŝ, is shown in Table 1. The full expressions for the LO helicity amplitudes are reported in Appendix A.1. 1 The operators studied here modify the partial width of the Z boson to quarks. The modification of the width is of order 1% if we saturate the c (1) ϕq , cϕu or c ϕd bounds, and smaller than 0.3% for c (3) ϕq . This modification can be reconciled with the SM branching-ratio predictions via (structural) cancellations with other dimension-6 operators that contribute to the same decays. These additional operators do not induce energy-growing corrections in the Zh process, but at most an overall rescaling. Since the sensitivity of our analysis comes from the shape of the differential distribution, rather than its normalization, the bounds we derive are not significantly affected by their exclusion. Table 1. High-energy behavior of the tree-level SM and BSM helicity amplitudes for the process qq → Zh.
For large partonic center-of-mass energy, √ŝ M Z , the SM amplitude with a longitudinally polarized Z boson is constant and dominates over the transversely polarized ones, which are suppressed by a relative factor of M Z / √ŝ . On the other hand, all the contributions corresponding to the dimension-6 operators (2.1)-(2.4) exhibit the same behavior. The amplitudes with a longitudinally polarized Z boson grow asŝ/Λ 2 , i.e., with the square of the center-of-mass energy, while, for a transversely polarized Z, the growth is only linear, of order √ŝ M Z /Λ 2 . A similar behavior characterizes the closely related W h process [15]. By contrast, in the W Z case, the transverse opposite-polarization channels are unsuppressed in the SM and constitute an important background for the BSM signal [2].
If one performs a fully inclusive analysis in the Z decay products, only amplitudes with the same Z polarization can interfere. In this case, the behavior of the squared SM amplitude and of the leading interference terms is 2 ϕq , M ϕu , M ϕd } and θ is the scattering angle of the Z boson. That is, all the leading BSM contributions can directly interfere with the leading SM amplitude. Thus, the main BSM effects are expected to be captured by exploiting the transverse-momentum distribution which is closely related to theŝ distribution.
However, stark differences in the size of the interference terms are present even though the behavior of the BSM amplitudes is the same for all effective operators (2.1)-(2.4). This can be seen by direct inspection of the analytic expressions (see Appendix A) and numerical analyses (see Tables 8 and 9). The only operator that leads to a sizeable interference is O ϕq , while all others suffer important suppressions.
For the right-handed operators, O ϕu and O ϕd , the interference terms are suppressed since they are proportional to the coupling of the Z boson to right-handed quarks, which is significantly smaller than the one to left-handed quarks. BSM effects linear in the Wilson coefficients are therefore suppressed with respect to the SM contributions and to the quadratic BSM terms, degrading the sensitivity. In fact we find that, for values of the   Wilson coefficients of order of the expected bounds (c ϕu ∼ c ϕd ∼ 2 × 10 −2 ), linear and quadratic BSM corrections are roughly of the same order.
For the weak isospin singlet operator, O ϕq , a (partially accidental) cancellation between the up-type and down-type quark contributions is present. This cancellation arises because the SM-BSM interference term is linear in the SM coupling of the quarks to the Z. Such coupling is proportional to T 3 − s 2 w Q, where T 3 is the weak isospin charge, Q the electric charge and s w is the sine of the weak mixing angle, and have opposite sign for up-type and down-type quarks. On the other hand, the leading BSM amplitudes (see Eq. (A.3)) coming from O (1) ϕq have the same sign for all quarks, leading to opposite-sign interference contributions for up-type and down-type quarks. 3 The suppression is further exacerbated by the relative weight of the valence-quark content of the proton, which makes the cancellation stronger. As a result, the sensitivity to O (1) ϕq is quite degraded and dominated by the terms that are quadratic in the Wilson coefficients.
A possible way to partially enhance the interference terms is to consider the differential distribution in the rapidity of the Zh system (or a correlated quantity, such as the Higgs rapidity). The rapidity distribution coming from up-type and down-type initiated processes is slightly different. In particular the uū parton luminosity is peaked at larger rapidity than the dd one (see for instance Fig. 3 of Ref. [28]). Separating the small rapidity and large rapidity regions, the cancellation can be partially undone. This is shown in Fig. 1 where the differential cross-section with respect to the Zh rapidity (normalized to the total SM cross-section) is plotted for different flavor assumptions. Note that only the interference term, which is linear in the Wilson coefficient, is shown. The solid, dash-dotted, double-dashed, and dashed curves correspond to the sum over the lightest 2, 3, 4, and 5 quarks respectively. This term changes sign around y Zh ∼ 2 because the relative contribution between the u-and d-quark initial states depends on y Zh as discussed above. The choice of N F depends on the flavor assumption imposed on the Wilson coefficient c (1) ϕq . The solid, N F = 2, curve sums over first-generation quarks which comprise the valence quarks. For any N F > 2, each down-type quark contributes positively in the central region (|y Zh | < 2) and each up-type quark contributes negatively. The size of the contribution progressively decreases the heavier the flavor. The resulting differences in the distributions could, in principle, be exploited to disentangle the flavor assumption given a clean decay channel with large statistics to sufficiently suppress the quadratic terms.
As we will discuss in Section 4, the impact of a rapidity binning on the analysis is small, mostly because the statistics in each bin is limited. Taking into account the rapidity distribution could however be more relevant for different final states with larger crosssection. For instance the Zh production process with the Higgs decaying into b quarks. We leave this analysis for future studies.
To conclude the discussion, we mention that interference between different Z-boson polarization channels can be recovered by considering differential distributions that also include the Z decay angles [29]. These effects, however, do not grow with energy and are small with respect to the main interference term which does. Taking into account the decay product distributions is also not useful to remove the suppression of the linear interference terms in c (1) ϕq , c ϕu and c ϕd . In fact all the interference amplitudes are controlled by the same combination of couplings and the same suppression patterns characterize all of them. We conclude that a sophisticated analysis taking into account the kinematics of the Z-boson decays can be avoided without a significant loss of sensitivity.

Additional contributions to the signal
So far we only considered the LO contributions to the Zh production channel. Additional contributions, with potentially different features, arise at the 1-loop level. In particular, the gg → Zh channel can show a very different dependence on the effective operators (2.1)-(2.4). Although new-physics contributions do not grow with the energy in this channel [19], a sizeable dependence on the effective operators might still be present. As shown in Ref. [30], the SM contribution to the gg → Zh channel is suppressed (especially at high energy) by a cancellation between the contributions from box-type and triangle-type diagrams. The presence of new-physics contributions can partially remove the cancellation, thus leading to sizeable deviations in theŝ distribution which grow at high energy. In fact we found that the corrections induced by O (1) ϕq for a value of the Wilson coefficient c (1) ϕq = 1.5 × 10 −2 (close to the bound we derive from our analysis) range from ∼ 4% in the low-energy bins up to ∼ 60% in the high-energy tail.
Extracting information from the gg → Zh channel can, however, be tricky for several reasons. First of all, its cross-section is much smaller than the leading qq → Zh one and, in the high-energy bins, where the corrections are stronger, it contributes to less than one event at the end of the FCC-hh operation (see Fig. 2). Special cuts would therefore be needed to enhance its visibility. Second, the gg → Zh process strongly depends on additional effective operators, whose determination could be not sufficiently precise to remove their effects. In particular, there is a strong dependence on deviations in the top Yukawa coupling. An uncertainty on its determination at the ∼ 1% level could easily overshadow the effects of the effective operators we are considering in this paper.
In the light of these difficulties, in our analysis, we will not exploit the new-physics dependence of the gg → Zh channel, ignoring its dependence on the Wilson coefficients and treating it as a background.
It is worth mentioning here another channel that can contribute to the signal. We characterize the final state with the Z decaying into neutrinos, by requiring a pair of photons and missing transverse momentum. A non-negligible fraction of events with this signature comes from the W h production channel with a leptonically-decaying W boson, in which the lepton is not detected. This happens in a significant fraction of the events in which the W decays into taus, while it is much rarer for decays into electrons and muons. The W h channel depends only on the O ϕq operator [2,15] with corrections that grow with the energy, thus enhancing the sensitivity to this operator.

Power counting considerations
We conclude this section with a short discussion about the expected size of the new-physics effects in generic BSM theories. All four operators considered here can be generated at tree-level via fermionic or bosonic weak isospin singlet and triplet states [31]. Therefore, from a power-counting point of view, no particular hierarchy exists between their Wilson coefficients. According to the SILH power counting [32], we find where g is the SM EW gauge coupling. This estimate is valid in theories in which new physics is weakly coupled or is not directly coupled to the SM fields.
If new physics is strongly coupled to the SM, the estimate becomes where g * is the typical size of the new-physics coupling.

Event generation and analysis
We simulated signal and background events with MadGraph5 aMC@NLO v.2.7.3 [33] using the NNPDF23 parton distribution functions [34]. Parton shower and Higgs decay into two photons were modeled using Pythia8.24 [35], while detector effects were estimated via Delphes v.3.4.1 [36][37][38][39][40][41] with its FCC-hh card. All the simulations were done using the SMEFTatNLO [42] UFO model [43], which implements the effective operators we consider in our analysis. Further details about the Monte Carlo generation can be found in Appendix B. The signal process is pp → Zh, with the Higgs decaying into a photon pair and the Z decaying either into a pair of charged leptons, Z → + − (with = e, µ), or into neutrinos, Z → νν. The channel in which the Z boson decays hadronically suffers from a much larger background and requires a more sophisticated analysis strategy that we leave for future work.

Z → νν channel
The main process contributing to this channel is qq → Z (→ νν) h (→ γγ). The signature is given by a pair of photons and sizeable missing transverse momentum. As we mentioned in Section 2.3, the additional process qq → W (→ ν) h (→ γγ) where the charged lepton is not reconstructed contributes to the same final state and must be considered as part of the signal, since it depends on the O ϕq operator. This has been found to be relatively unlikely for a W to decay into electrons and muons, but happens in a non-negligible fraction of events with W → τ ν τ (see Section 3.3). We simulated these channels at LO in QCD and we took into account NLO QCD and EW corrections through k-factors. NLO corrections turn out to be sizeable, but with a weak dependence on new-physics. It is thus a good approximation to rescale the cross-section in each bin by the SM NLO k-factors which are given in Table 6. We also take into account the gg → Z(→ νν)h(→ γγ) production channel at LO but neglect new-physics contributions as discussed above.
The main backgrounds for this channel come from the processes Z (→ νν) γγ and W (→ ν) γγ, where the charged lepton is missed and the photon pair is non-resonantly produced. Additional contributions from events in which one or two jets fake a photon are negligible if we assume a P j→γ = 10 −3 jet-to-photon fake rate [44,45] (see also discussion in Appendix B of Ref. [15]). Consequently, we do not consider processes with less than 2 photons as part of the background. This is valid for both Z decay channels. We also checked that the jγγ process in which missing transverse momentum comes from showering and detector effects, constitutes a huge background at low energy, but becomes negligible when a 200 GeV cut on the missing transverse momentum is imposed.
The backgrounds were generated at NLO in QCD. This is particularly important for the W (→ ν) γγ and Z (→ νν) γγ channels. In the former, the cross-section is fully dominated by the NLO contribution with a real emission (which is a factor ∼ 17 larger than the LO contribution). Similarly, the NLO contributions for the Z (→ νν) γγ channel are a factor ∼ 2 larger than the LO ones. The enhancement is partly due to the fact that at NLO an additional channel with an initial gluon instead of a sea quark opens up. As we will see in Section 3.3, the presence of a large component of events with an additional jet significantly modifies the kinematic configuration of the events, providing efficient handles to distinguish the background processes from the signal.
The main signal process in this channel is qq The main background process is Z(→ + − )γγ, with a small additional contribution from processes in which the lepton pair arises from the splitting of a virtual photon. This contribution can be suppressed by restricting the di-lepton invariant mass to be around the Z-pole. As for the neutrino channel, we generated the signal at LO QCD, applying suitable k-factors to take into account NLO QCD and EW effects (see Table 6), while the background was simulated at NLO in QCD. Clearly, the NLO corrections follow a pattern completely analogous to the one we discussed in the neutrino channel. We also took into account at LO the gg → Z(→ + − )h(→ γγ) channel, considering it as a constant background.
Finally, we mention that the decay of the Z boson into a pair of taus, when both taus decay fully leptonically, can also contribute to the + − channel. We estimate that this process would increase the signal and background cross-sections in a similar way, leading to at most a 5% increase. The impact of such change on our results would be negligible, therefore we do not include this channel in our simulations.

Selection cuts
In order to reconstruct the signal we require the presence of at least two photons with p T > 50 GeV and |η| < 6. In addition, at least one pair of photons must have an invariant mass in the range 120 GeV < m γγ < 130 GeV. If more than one pair of photons fulfill that requirement, we choose the pair with the smallest angular separation, ∆R γγ = (∆η γγ ) 2 + (∆φ γγ ) 2 . 4 In the Z → νν channel, we veto events with muons or electrons in the kinematic region defined by p T > 30 GeV and |η| < 6. In the Z → + − channel, we only accept events which contain 2 leptons of opposite charge, each with p T > 30 GeV and |η| < 6. The invariant mass of the lepton pair m l + l − is required to be in the range [81, 101] GeV to ensure the reconstruction of the Z boson mass.
Finally, we impose a maximum cut on the transverse momentum of the reconstructed Zh system, p Zh T < p Zh T,max . This cut is motivated by the fact that a large transverse momentum is typically the signal of a recoil against a hard QCD jet. Such events are more likely to come from background processes than from the signal. A summary of the selection cuts and of the values of p Zh T,max is given in Table 2. To conclude this discussion, let us mention that we also checked the usefulness of an upper cut on the angular separation of the photon pair, ∆R γγ , and of the lepton pair, ∆R + − . These cuts are motivated by the fact that the di-photon and di-lepton pairs, coming from boosted objects, typically have smaller ∆R than the non-resonant background processes. We however found that, after the cuts on the invariant masses, the ∆R distributions for signal and backgrounds were very similar, so that a cut on the angular separation does not improve the fit.

Binning
In our analysis, we consider the double-differential distribution in the center-of-mass partonic-energy and in the rapidity of the events (see Section 2). As can be seen from Eq. (2.6), the signal events are mostly emitted in the central scattering region. We can therefore approximately trade the event energy for the transverse momentum of the Z or Higgs. Both of these can be reconstructed in the two decay channels we consider. We choose to bin in the minimum p T of the two bosons where, in the neutrino channel, p Z T is identified with the missing transverse momentum. This choice is useful to select hard events in the Zh center-of-mass frame and, simultaneously, to remove the jγγ background where the missing energy comes from soft radiation. We use five bins in p T,min , whose boundaries are given by p T,min ∈ {200, 400, 600, 800, 1000, ∞} GeV . (3.2) In the Z → + − channel, we also use a simple binning in the rapidity of the Zh system y Zh , namely |y Zh | ∈ [0, 2], [2,6] , while in the Z → νν channel, since the rapidity of the Zh system cannot be determined, we bin in the rapidity of the Higgs y h instead, which is strongly correlated with y Zh . The bin definitions depend on p T,min as follows,

Cut efficiencies
We show in Tables 3 and 4 the cutflow analysis for the different signal and backgorund processes in the Z → νν and Z → + − channels. In order to focus on events that have a good sensitivity to new-physics, in obtaining these results we only considered high-energy events that satisfy the parton-level cut p Z,W T > 400 GeV. Notice that the initial phase space of the different processes does not exactly coincide, due to small differences in the generation-level cuts (see Appendix B and Table 7). The results presented in the tables are however still useful to understand the efficiency of each cut in suppressing the backgrounds.
As expected, the cut on the invariant mass of the photon pair, m γγ , is very efficient in reducing all backgrounds with non-resonant photon pairs, i.e., W γγ and Zγγ. Analogously, the cut on the lepton-pair invariant mass m + − is helpful in reducing the + − γγ  Table 4. Cutflow for the processes in the Z → + − channel. The acceptance region for charged leptons and photons is defined in the text.
background component in which the lepton pair comes from an off-shell photon line. Such channel becomes completely negligible after the cuts, so that we did not include it in Table 4, where only the Zγγ → + − γγ channel is listed. It is interesting to notice that, in the neutrino channel, the charged-lepton veto reduces the cross-section of the W h and W γγ channels by only 60 − 70%. This behavior, which might seem surprising at first sight, is mainly due to the W → τ ν τ decay channel, in which the τ decays hadronically. On the other hand, when the W decays to light leptons, the cut is much more efficient, with a reduction of the events of order 96% for decays with an electron and 98% for decays with a muon.
A large fraction of the events with hadronically decaying taus is however removed by the cuts on p T,min and p Zh T . The reason for this is twofold. First, the jets from the τ decay carry a sizeable fraction of the τ momentum, thus creating an unbalance between the missing transverse momentum and the photon pair p T . And, second, in the W γγ channel, a vast majority of the events come from the NLO contribution with an extra hard jet. In such events a large transverse momentum p Zh T is present. The p T,min cut is also quite efficient in reducing the Zγγ backgrounds both in the neutrino and lepton channels. This, again, is due to the fact that many events (∼ 65%) come from the NLO contribution with an additional hard jet.  Tables 8 and 9, together with the dependence of the signal cross-section on the Wilson coefficients.
In the neutrino channel, it is remarkable that the W h channel gives a significant contribution to the signal with a number of events of the same order of magnitude as in the Zh channel. This is true in spite of the fact that the W h channel has a much lower cut efficiency with respect to Zh, as can be seen from Table 3. The suppression is however partially compensated by the significantly higher initial cross-section. Moreover the dependence on c  Tables 8 and 9.
The main background in the neutrino channel is Zγγ → ννγγ, followed by W γγ. In most of the p T,min bins, the backgrounds are much smaller than the signal, making the analysis almost background free. The only exceptions are the first bin, where however the dependence on new physics is small, and the last one, where the number of events is quite limited. Finally the gg → Zh channel is also quite small, so that its dependence on new physics can be safely neglected.
In the charged lepton channel, the Zγγ → + − γγ background is much smaller than the signal in all bins. Analogously to the neutrino channel, the gg → Zh has an almost negligible impact, especially in the higher p T,min bins.
Before discussing the results, we make a comment on the impact of the flavor-universality assumption for the new-physics contributions. The main signal channels Zh and W h are dominated by the production through initial states containing first-generation quarks. Contributions from second-generation quarks and from the bottom are suppressed by the proton parton distribution functions, and together give ∼ 20% of the whole cross-section. 5 Thus, the analysis we perform approximately captures the case in which new physics couples only the first generation, or the case of U(2) flavor symmetry in the first two generations.
The situation is potentially different for the gg → Zh channel. In this case new-physics corrections to the top-quark couplings can strongly affect the dominant loop diagrams. In the flavor-universal scenario, the corrections to the top couplings are related to the ones to the light generations, thus they are relatively small. We checked that, for values of the Wilson coefficients of the order of the bounds we find (|c ϕq |, |c ϕu |, |c ϕd | ∼ 2 × 10 −2 ), the corrections to the gg → Zh cross-section are below 5% in the first three bins and become ∼ 30% in the fourth bin and ∼ 60% in the last bin. Given the small number of events in this channel, the new-physics impact on the fit is clearly negligible.
The situation could change in the case in which the flavor-universality assumption is relaxed and the top-coupling modifications are allowed to be larger. In this case a dedicated analysis strategy exploiting other processes sensitive to the top couplings, in addition to the gg → Zh channel, would be required. We notice, however, that relaxing flavor universality could easily induce large flavor-violating effects, unless flavor symmetries or alignment assumptions are imposed.

Results
In this section, we present our projection of the bounds on the Wilson coefficients. We first focus, in Section 4.1, on the O ϕq operator than on the others, as explained in Section 2.2, so that a much more stringent bound is expected. Moreover, in many new-physics scenarios we expect the size of the various operators to be comparable (see Section 2.4), in which case the constraints on O   To give an idea of the impact of systematic uncertainties, we considered three benchmark scenarios, with 1%, 5%, and 10% uncorrelated systematic error in each bin. We notice that the bounds for positive and negative values of the Wilson coefficient are quite similar. As can be checked from the results in Tables 8 and 9, this is due to the fact that the new-physics contributions are dominated by the linear interference terms with the SM. 5 The relative contribution of each qq initial state to the total cross-section at LO in the bin 400 GeV < p T,h < 600 GeV is: u(42%), d(38%), s(9.2%), c(5.7%), b(4.7%).
Interestingly, the bounds in Eq. (4.1) are competitive with (and actually slightly better than) the ones obtained with a similar one-operator fit in the W h → νγγ production channel [15], which for 5% systematic error gives c  For comparison, the projections obtained for the leptonic W Z channel at FCC-hh, assuming 5% systematics, give a bound [2] FCC-hh (20 ab −1 ) c which is only slightly more stringent than the bound we find combining the W h and Zh channels. Further combination with the W Z bounds can therefore lead to a significant improvement in the sensitivity. Other current and projected bounds at 95% C.L. on c  We also quote in parentheses the bound from one-operator fits. All these bounds are much weaker than the ones from the W h, Zh and W Z channels, with the exception of the FCC-ee and CEPC one-operator fits.

Full analysis
The bounds from the four-operator fit, i.e, with O ϕq , O ϕu and O ϕd , profiled over three operators at a time are reported in the column labeled "Profiled Fit" in Table 5. These results clearly show that, even in a full fit, the sensitivity to the O (3) ϕq operator is much higher than to the other operators.
The row labelled c Table 5 corresponds to the combination with the bounds from the one-operator fit in the W h channel presented in Ref. [15]. Notice that, in the W h channel, the operator O (3) ϕq is the only one that grows quadratically with the energy, which justifies our choice of the one-operator fit. This combination leads to an improvement by a factor ∼ 2 in the bounds on c     ϕq (+W h [15]) provides the bounds obtained from the combination of our analysis with the one for the W h channel presented in Ref. [15]. Left column: Bounds from the global fit, profiled over the other coefficients. Right column: Bounds from a one-operator fit (i.e. setting the other coefficients to zero). the fact that the bounds from single-operator fits on c (1) ϕq , c ϕu and c ϕd are nearly equal to the ones coming from the global fit (see "One-Operator Fit" column in Table 5). On the contrary, the determination of c ϕq . In Fig. 3, we provide the fits in the (c ϕq terms) can be seen, which however disappears in the profiled fit. In the right panel of the figure we also show how the fit is modified by combining with the W h analysis of Ref. [15]. The main impact of the combination is an improvement of the bounds on c  Finally we can compare the bounds we obtained with the expected sensitivity at the LHC and other proposed future colliders (our fits always implicitly also rely on LEP measurements that are crucial in particular to size the SM input parameters). The current and ϕq , c ϕu and c ϕd from one-operator fits as functions of the maximalinvariant-mass cut M . The bounds correspond to 95% C.L. (∆χ 2 = 3.84). The dashed, solid and dotted lines show the bounds for 1%, 5% and 10% systematic errors. The current bounds and the projections for some future hadron colliders are also shown. For c (3) ϕq , we show the LHC Run 1 bound from Ref. [49] and the projections from the W Z channel at HL-LHC from Ref. [2]. For the rest of operators, we show the global LHC data fit bound from Ref. [48] and the 1-operator fit at HL-LHC from Ref. [47]. LHC bounds are shown in orange and HL-LHC ones are shown in dark green. The dashed gray and solid black lines show the values of the Wilson coefficient expected in weakly-coupled (c ∼ g 2 /(4M 2 )) and strongly-coupled (c ∼ (2π) 2 /(M 2 )) new physics models [2]. We see that, for all operators, our analysis provides bounds that are competitive with the ones expected from global fits at other future colliders. On the other hand, if oneoperator fits are considered, FCC-ee and CEPC will surpass our bounds on all four operators by roughly one order of magnitude.
The bounds from single-operator fits as functions of the maximal invariant mass of the events used in the analysis are shown in Fig. 4. This kind of analysis is very useful in testing the validity of the EFT description and understanding what kind of theories can be tested. We see that, for all the operators, the bounds saturate for M ∼ 5 TeV, signalling that our bounds are valid for a cut-off Λ 5 TeV. One can also appreciate the fact that only the bounds on the O ϕq operator are strong enough to test the region of weakly-coupled new physics, whereas the other operators can mostly test theories in which the new dynamics is strongly coupled.

Connection to aTGCs and Universal Theories
So far we used our analysis to set bounds on the O ϕq , O ϕu and O ϕd operators can be rewritten in terms of vertex corrections, δg Zq L,R , and aTGC, δg 1z and δκ γ by adopting the Higgs basis [50,51]: where c w , s w and t w are the cosine, sine and tangent of the weak mixing angle respectively. From Eq. (4.8) we see that if the δg Zq L,R coefficients are negligible, one can recast our diboson bounds on c ϕq , c ϕu , c ϕd to constraints on δg 1z and δκ γ . This is particularly justified for universal theories where δg Zq L,R correspond to combinations of the oblique parameters S, T, W and Y [2,6], which are expected to be highly constrained through various measurements at FCC-ee and FCC-hh. With these assumptions, the combined analysis of the Zh and W h channels gives the following constraints (for 5% systematic uncertainty) where the bounds in parenthesis are obtained through one-operator fits. Another way to estimate the impact of diboson searches at FCC-hh is to compare their reach with the sensitivity at future lepton colliders. As an illustrative example we show, in Fig. 5, the expected bounds on the aTGC parameters δg 1z and δκ γ for CEPC and FCC-ee. These bounds are obtained through a global fit, which takes into account 18 effective operators in flavor-universal theories [47]. 6 We also show in the same plot the impact of the combination of the CEPC and FCC-ee constraints with the FCC-hh fit of diboson channels. For the FCC-hh fit we combine the projections for W Z → ν + − [2], W h → νγγ [15], and the two channels Zh → + − γγ and Zh → ννγγ (this work), and we assume 5% systematic uncertainty.
We find that, once the diboson channels at FCC-hh are included, the projected bounds on δg 1z improve by a factor ∼ 2 (3) with respect to the bounds at FCC-ee (CEPC). On the other hand the bounds on δκ γ are only marginally affected. This behavior is due to Figure 5. Projected 95% C.L. bounds on the anomalous Triple Gauge Couplings δκ γ , δg 1z for flavor universal theories at future colliders. In dark and light orange, we show the expected bounds at FCC-ee and CEPC [47]. In dark (light) blue, we show the combination of the FCC-ee (CEPC) fit with the projections at FCC-hh for the diboson channels W Z → ν + − [2], W h → νγγ [15], and the combined Zh → + − γγ and Zh → ννγγ (this work). In the FCC-hh projections we assume 5% systematic uncertainty. the fact that diboson channels have a very good sensitivity to c (3) ϕq , which only depends on δg 1z (see Eq. (4.8)), while information on δκ γ can only be extracted from c (1) ϕq , c ϕu and c ϕd , whose expected bounds are much weaker.
For completeness we report in the following the allowed 95% C.L. regions for δg 1z ,

Summary and conclusions
The next generation of hadron colliders, thanks to the enhanced cross-sections and the high integrated luminosity, will allow us to revisit several electroweak processes in a cleaner environment, bringing advantages for precision measurements. An interesting example is V h production, which can only be studied at HL-LHC through the h → bb decay, but becomes accessible at FCC-hh also in the very clean channel h → γγ. In a previous paper [15], we exploited the latter channel in W h production, showing how it can be used to test deviations in the W -boson couplings to quarks with high accuracy. In the present work, we extended the analysis to a closely related channel, Zh production, considering the h → γγ decay channel together with a Z-boson decay into a pair of charged leptons or neutrinos. In spite of its smaller cross-section, the Zh channel is extremely interesting since it can be used to probe a larger set of new-physics effects, including deviations in the Z couplings to quarks that cannot be tested in the W h channel. In the SMEFT framework, and focusing on dimension-6 operators that induce aŝ growth in the amplitude, the Zh production process is sensitive to four "primary" operators [2], which in the Warsaw basis correspond to O ϕq , O ϕu , and O ϕd . Binning in transverse momentum and combining both aforementioned Z decay channels allowed us to obtain a sensitivity on O (3) ϕq that is competitive with other processes that have higher cross-sections, like W h. Our estimates show that the bounds at FCC-hh can significantly surpass the precision achievable at HL-LHC and at FCC-ee. Furthermore, the combination of the Zh and W h channels can improve the bound by roughly 20% and further improvement could be achieved via a combination with the W Z channel. This shows the importance of a comprehensive study of all processes available at future colliders in order to correctly assess their potential.
The sensitivity to the O ϕq operator turns out to be significantly smaller, due to a (partially accidental) cancellation between the contributions from up-type and down-type quarks. We showed that this cancellation can be partially overcome by implementing a second binning in the rapidity of the Zh system (or the rapidity of the Higgs boson when the former cannot be computed). This double binning exploits the differences in the rapidity distribution due to the parton distribution functions of the various generation-level quarks. Although the improvement is limited in the final state we considered in this paper, this strategy could be useful for final states with higher cross-section (for instance the h → bb channel) or for analogous processes. We also stress that the pattern of deviations in the rapidity distribution depends on the flavor structure of the new-physics effects (see Fig. 1), thus it could potentially be a way to disentangle different flavor hypotheses. The sensitivity we get to the O (1) ϕq operator is one order of magnitude better than the one at HL-LHC and is competitive with the one achievable at future lepton colliders (CLIC/ILC, CEPC and FCC-ee).
Finally, our sensitivity to O ϕu and O ϕd is limited by inherent characteristics of the operators, whose contributions only interfere with the SM amplitudes due to the relatively small couplings of the right-handed quarks with the Z boson. The expected bounds are more than one order of magnitude better than the ones at HL-LHC and competitive with ϕq , c ϕu and c ϕd . In blue, our combined bounds from Zh → (νν/ + − ) γγ and W h → νγγ at FCC-hh with 30 ab −1 for different systematics and computed from a four operator fit. In all cases, the black lines with a triangle on top represent the bound from a one-operator fit instead. In light yellow, the current LEP [46] bound for c (3) ϕq . In light green for c (3) ϕq , the run-1 LHC [49] bounds. In medium green, the current bound on all the operators from a global fit [48]. In dark green, the projections from a global fit at HL-LHC [26,47]. In light, medium and dark orange, the projected bounds on the operators from a global fit at CLIC, CEPC and FCC-ee respectively [47]. FEPC stands for Future Electron-Positron Colliders. the ones from global fits at future lepton colliders. Regarding the bounds on O ϕu and O ϕd , we note that, due to the suppression of the interference with the SM amplitudes, the constraints are mostly driven by the square of the BSM contributions. This might cause some limitation in their interpretability within the EFT formalism.
A summary of the projected 95% C.L. bounds on the four operators we considered is shown in Fig. 6. The blue bars correspond to the constraints derived from the profiling of a four-operator fit. On the other hand, the horizontal bars with a triangle indicate the bound obtained from a fit including one operator at a time. In both fits, we considered three possible values for the systematic uncertainties: 1% (lighter shading), 5% (medium shading), and 10% (darker shading). The systematic uncertainty has a sizeable effect only on the bound for O (3) ϕq . The 5% scenario is comparable to the present LHC systematics for similar processes, therefore it could be considered as a conservative estimate, while the 10% benchmark is most probably an over-pessimistic one.
Many directions in the assessment of the precision-measurement potential of future hadron colliders could still be explored. Regarding the Higgs-associated production channels (W h and Zh), an interesting direction to follow is the study of the hadronic decay channels, in particular the h → bb decay which can offer a boost in cross-section, however at the price of significantly larger backgrounds. Hadronic decays of the W and Z bosons could also be considered. In the broader context of diboson channels, the W Z and W W production processes can also be exploited to test a similar set of new-physics effects. The former channel has been so far investigated mainly in the fully leptonic final state, while other decays could also be worth considering. W W production is instead much more challenging and very few studies are already available. Finally, we stress that such analyses would make a global and combined fit of all the diboson channels possible, and we can expect this to further improve the sensitivity to new physics fully exploiting the potential of FCC. It is worth noting that a combination of all channels could also be beneficial in reducing systematic and theory uncertainties, which are partially correlated across different processes.  Table 6. NLO k-factors for the main signal processes. Each entry shows separately the QCD and QED contributions to the k-factors. The accuracy on the determination of the k-factors is of order few × 10 −2 .
A.2 Squared amplitudes and interference terms for qq → Zh The squared SM amplitude and the SM-BSM interference terms are given separately below. For convenience, we define a function that depends on the scattering angle and is common among all squared and interference terms,  The k-factors in the various p T,min bins are listed in Table 6 in the format 1 + (k qcd − 1) + (k qed − 1). As one can see, they give an enhancement of the cross-section of up to 50%. QCD corrections, which enhance the cross-section, typically dominate. On the other hand, QED effects, which tend to lower the cross-section, are subleading in the low-energy bins, whereas they become comparable to the QCD ones at high energy.
The background processes were simulated at NLO in QCD, with the exception of gg → Zh, which was simulated at LO. QED corrections have not been applied to the backgrounds. This is a conservative choice since they reduce the cross-sections.
Generation cuts We applied a set of cuts at generation level to increase the number of Monte Carlo events after the selection cuts. These generation cuts are shown in Table 7, where the second (third) column shows the values used for all the processes relevant for the Z → νν (Z → + − ) channel. Furthermore, we generated the events in 6 exclusive bins in the p T of the gauge boson. We notice that, due to initial state radiation effects, a sizeable migration of events between generation and p T min bins is present. Tables 8 and 9, we show the fits of the signal and background cross-section in the various bins as a function of the c   Table 9. Number of expected signal events as a function of the Wilson coefficients (with Λ = 1 TeV) and background events in the Z → + − channel at FCC-hh with 30 ab −1 integrated luminosity. The Monte Carlo errors on the fitted coefficients, when not explicitly specified, are at most of order few %.