Hadron Collider Probes of the Quartic Couplings of Gluons to the Photon and $Z$ Boson

We explore the experimental sensitivities of measuring the $gg \rightarrow Z \gamma$ process at the LHC to the dimension-8 quartic couplings of gluon pairs to the $Z$ boson and photon, in addition to comparing them with the analogous sensitivities in the $gg \to \gamma \gamma$ process. These processes can both receive contributions from 4 different CP-conserving dimension-8 operators with distinct Lorentz structures that contain a pair of gluon field strengths, $\hat G^a_{\mu \nu}$, and a pair of electroweak SU(2) gauge field strengths, $W^i_{\mu \nu}$, as well as 4 similar operators containing a pair of $\hat G^a_{\mu \nu}$ and a pair of U(1) gauge field strengths, $B_{\mu \nu}$. We calculate the scattering angular distributions for $gg \rightarrow Z \gamma$ and the $Z \to \bar f f$ decay angular distributions for these 4 Lorentz structures, as well as the Standard Model background. We analyze the sensitivity of ATLAS measurements of the $Z(\to \ell^+\ell^-, \bar \nu \nu, \bar q q)\gamma$ final states with integrated luminosities up to 139 fb$^{-1}$ at $\sqrt{s} = 13$ TeV, showing that they exclude values $\lesssim 2$ TeV for the dimension-8 operator scales, and compare the $Z \gamma$ sensitivity with that of an ATLAS measurement of the $\gamma \gamma$ final state. We present combined $Z \gamma$ and $\gamma \gamma$ constraints on the scales of dimension-8 SMEFT operators and $\gamma \gamma$ constraints on the nonlinearity scale of the Born-Infeld extension of the Standard Model. We also estimate the sensitivities to dimension-8 operators of experiments at possible future proton-proton colliders with centre-of-mass energies of 25, 50 and 100 TeV, and discuss possible measurements of the $Z$ spin and angular correlations.


Introduction
In the Standard Model (SM) of particle physics, the triple and quartic gauge couplings are fixed by gauge symmetry. Measuring these couplings can test not only whether the gauge symmetry is realized linearly or nonlinearly in the low-energy effective theory of the electroweak symmetry-breaking sector [1] but also provide an interesting way of looking for possible new physics beyond the SM [2,3]. Hence the search for anomalous gauge couplings is one of the priority measurements for LHC and possible future colliders. Many studies have been made of the present and prospective experimental sensitivities to triple gauge couplings (TGC) and quartic gauge couplings (QGC) between electroweak SU(2) L ×U(1) Y gauge bosons. Quartic couplings between these and the SU(3) c gluons are absent in the Standard Model (SM), but are also allowed in the Standard Model Effective Field Theory (SMEFT) at the level of dimension-8 operators. These have not been studied to the same extent as dimension-6 SMEFT operators (see [4,5] and references therein), though there have recently been studies of dimension-8 operators that generate QGCs between photons [6] and between photons and gluons [7], as well as those that generate neutral TGCs [8,9].
The dimension-8 operators that generate QGCs also arise from loop diagrams in the SM, via extensions of the original calculations by Heisenberg and Euler [10,11], and may be generated by the exchanges of massive axion, dilaton or spin-2 resonances. One particular combination of dimension-8 interactions arises in Born-Infeld theory [12][13][14][15]: where the index λ runs over the 12 generators of the SM SU(3) c ×SU(2) L ×U(1) Y gauge groups. The parameter β ≡ M 2 is the Born-Infeld nonlinearity scale associated with high-scale physics. The Born-Infeld extension of SM (BISM) may have deep roots in M-theory-inspired models, where it is related to the separation between branes [14]. Since all the SM gauge group factors appear in Eq. (1.1), it generates quartic couplings of gluons to electroweak gauge bosons (gluonic QGCs, gQGCs) [7].
The cleanest strategy for searching for the photonic and gluonic QGC operators is to study processes that do not receive SM or dimension-6 SMEFT contributions. For example, the diphoton final state generated by light-by-light scattering, γγ → γγ, provides a very clean probe of the photonic QGCs [6]. A first LHC constraint on light-by-light scattering was provided by an ATLAS measurement in heavyion collisions [16] (see also [17]). Its rate was found to be consistent with the Heisenberg-Euler prediction [10,11], allowing lower limits O(100) GeV to be set on the scale of a Born-Infeld extension of QED and other possible dimension-8 SMEFT interactions [6]. Recently, the CMS and TOTEM collaborations updated the lower limits on the basis of a search for the exclusive production of highmass γγ final states in pp collisions at the LHC [18]. Diphoton final states generated by gluon-gluon scattering, gg → γγ, provide clean probes of gluonic QGCs [7], and the 13 TeV data of ATLAS with 37 fb −1 [19] enabled lower limits 1 TeV to be placed on the scales of the gQGC operators [7].
In this paper we study in more detail the present and prospective future experimental sensitivities to the gQGC operators involving gluons and pairs of neutral electroweak bosons, Z and γ. We present a first analysis of gg → Zγ scattering and compare this with its gg → γγ counterpart. We analyze ATLAS measurements of Zγ production followed by Z → + − [20],νν [21] andqq [22] decays, with up to 139 fb −1 of luminosity at √ s = 13 TeV. The scale of the dimension-8 gQGC operators can be constrained up to 2 TeV, higher than obtained in the updated analysis of the diphoton channel that we also present in the current paper. However, we also show that whereas measurements of the Zγ final state do not constrain the BISM scale, the γγ process constrains its scale to 5 TeV. Interestingly, we note that a combination of measurements of the Zγ and γγ final states could in principle disentangle the contribution of different operators. Finally we display the combined sensitivities of the Zγ and γγ channels, and show how the sensitivities to dimension-8 operators could be increased in the future by measurements at possible proton-proton colliders with centre-of-mass energies of 25, 50 and 100 TeV.
The structure of the paper is as follows. Sec. 2 reviews the dimension-8 gQGC operators, especially those relevant to the associated production of Z and γ, as well as γγ pair-production. Then Sec. 3 studies the kinematics and differential cross sections of the gg → Zγ process at the LHC. In particular, we show how the unitary constraint can be implemented consistently. Following these preparations, we investigate the experimental constraints on the gQGC operators at the LHC in Sec. 4 and future colliders in Sec. 5. Then, in Sec. 6 we present the helicity amplitudes for the partonic process gg → Zγ and analyze the corresponding angular distributions as well as spin correlations in Z boson decay, and discuss the possible interest of measurements of the Z spin and angular correlations. Our conclusions can be found in Sec. 7.

Dimension-8 Gluonic QGC Operators
Although there are many possibilities for new physics at higher scales above the EW one, the low energy effective field theory (EFT) should be subject to the SM SU(3) c ×SU(2) L ×U(1) Y gauge symmetries. A convenient way to take these symmetries into account is the SMEFT [23], which includes systematically all the allowed interactions with mass dimension d > 4. The extra dimensions are compensated by inverse powers of a mass scale M that is associated with heavy new physics. The gQGC operators appear at dimension-8 level with 1/M 4 suppression [7], In order to respect the SM gauge symmetries, the gauge bosons should appear via their field strengths, such as G a µν for gluons together with W i µν and B µν for EW gauge bosons. Since gluons carry QCD color, denoted by the a superscript of G a µν , the gQGC operators must contain an even number of gluon field strengths, such as G a µν G a,αβ , so as to be colorless. The same thing applies for the SU(2) L ×U(1) Y gauge boson field strengths, for example W i µν W i,αβ and B µν B αβ . Another symmetry to be imposed is Lorentz invariance. There are four different Lorentz-invariant contractions, as shown above. Hence we must consider eight gQGC operators in total.
The total and differential cross sections are largely determined by the Lorentz structure of the gQGC operators. The eight operators can be classified into four pairs {O gT,(0,4) , O gT, (1,5) , O gT, (2,6) , (3,7) }, each with a same Lorentz structure. The W i and B fields do not correspond to the physical mass eigenstates, which are the photon A and Z boson fields, W 3 = c w Z + s w A and B = c w A − s w Z where (c w , s w ) ≡ (cos θ W , sin θ W ) are the cosine and sine functions of the weak mixing angle θ W . Each pair W i W i and BB contains AA, AZ and ZZ contributions with coefficients determined by θ W : We see that their contributions to gg → Zγ scattering via the ZA combination differ only by a sign. Consequently, they produce exactly the same event distributions, and therefore yield identical cross sections for common values of the dimension-8 cut-off scale.
Only 4 independent operators contribute to gg → Zγ. On the other hand, in the case of gg → γγ, one needs to study separately the 8 operators, as done in [7], and the contributions of BB operators are larger than those of their i W i W i counterparts by a factor cot 4 θ W . Thus the γγ and Zγ channels in gluon-gluon scattering are complementary, and their combination can help to disentangle the W i W i and BB components within each pair.
For comparison, the Born-Infeld extension of the SM shown in Eq. (1.1) contains only one linear combination of the eight gQGC operators, where the traces are over the Lorentz and SM gauge group indices. The first term involves only CPeven pairs GG and W W of gluon and EW field strengths while the second is composed of the CP-odd combinations G G and W W . The combinations of CP-odd operators can be rewritten in the form of Eq. (2.1), as follows: Tr(G G)Tr(W W ) = 4Tr(GW GW ) − 2Tr(GW )Tr(GW ), (2.4) and similarly for Tr(G G)Tr(B B). Hence the BISM contribution to gQGC is simply a linear combination of six gQGC operators in the basis considered above, namely Most importantly, these six operators share a common cut-off scale M and their coefficients are fully correlated.
It is important to observe that, since the different gauge sectors appear with equal coefficients in the Born-Infeld lagrangian Eq. (1.1), the dimension-8 couplings of the EW gauge bosons to gluons take the diagonal form i W i W i + BB. Hence, because of the opposite signs and equal magnitudes of the ZA terms in Eq. (2.2), the BI extension of the SM does not contribute to Zγ production, but only to γγ and ZZ production. As we shall see in Sec. 4, the sensitivity to the scale of the BI extension of the SM at the LHC via gg → γγ scattering is considerably greater than that obtained from the study of γγ → γγ scattering.

Properties of the gg → Zγ and γγ Cross Sections
As discussed in Sec. 2, dimension-8 contributions to gluon-gluon collisions can yield γγ and Zγ final states. Both can provide a clean signal for probing gQGC operators, whilst the BISM can be probed only by the γγ channel. This section summarizes the properties of the total and differential gg → Zγ and γγ cross sections, discussing the basic features and laying the basis for the analysis of the various Z decay modes in Sec. 4.

Total Cross Section and Energy Dependence
When regarded as a probe of new physics beyond the SM, the sensitivity of the search for a gQGC signal is limited by the possible event rates. A conservative search strategy starts by using the total cross section. Since each gauge field strength (G a µν , W i µν , and B µν ) contains one derivative and hence contributes one factor of momentum p to the scattering process, the four field strengths in a dimension-8 operator yield a scattering amplitude M with a p 4 dependence. Taking into account the flux prefactor that is ∝ 1/ŝ for high-energy scattering, we find that the cross sections given by the gQGC operators grow asŝ 3 ∼ p 6 . Hence the signal cross sections generated by the O gT,i operators are rapidly-rising functions of the center-of-mass energyŝ of the gluon pair: where x Z ≡ m 2 Z /ŝ. This effect enhances the signal event rate significantly at higher energy, as shown in Fig. 1, where we see that the signal grows very fast whereas the SM background decreases ∝ 1/ √ŝ . Thus one expects to find greater sensitivities at higher-energy colliders, as we discuss later. The steadily growing cross sections finally violate the unitarity bound at some center-of-mass energy scale [7]. This problem can be addressed by introducing momentum-dependent couplings or unitarization procedures [24], or by estimating the maximum center-of-mass energy allowed by unitarity [25]. Following the procedure adopted in [6] and [7], here we respect the unitarity constraint by assuming that the cross section falls asŝ −1 once the "tree unitarity" requirement [26] is satisfied, after which the energy scaling is the same as that of the major irreducible background due toqq → γγ. We find that this unitarity bound σ ∼ 1/m 2 Zγ is saturated at the gluon-gluon center-of-mass energies for i = (0, 4), (1,5), (2,6), (3,7), in the limit x Z → 0. Forŝ = 1 TeV 2 , the ratio x Z 0.01 is generally negligible.
The sensitivity to the scale of new physics is quite stable with respect to the choice of unitarization scheme. For example, an aggressive unitarization scheme that cuts off the signal events above the saturation point would reduce the total event rate by a factor of about 2. But the dependence on the new physics scale M i is eighth power, σ ∝ 1/M 8 i , so the sensitivity is reduced by only 10% at most. In following discussions, the same unitarization scheme σ ∼ 1/ŝ is applied for both the current LHC-13 TeV with real data and the future experiments with projected pseudo-data.
The total cross sections for gg → Zγ and γγ are summarized in Fig. 1 as functions of √ŝ , with M i = 1 TeV for illustration. The signal cross sections are much larger than the SM background for O gT, (1,5) O gT, (2,6) O gT, (3,7) 0.
1000    (1,5), (2,6), (3,7) for the gg → Zγ process; √ŝ i = (4.71, 6.21, 6.51, 6.88, 3.49, 4.60, 4.82, 5.10)M i for i = 0, . . . , 7 for the gg → γγ process; and √ŝ i = 3.26M i for the BI model. Even when the unitarity bound is imposed for gg → Zγ, replacing its rapidly-growing cross section with the SM 1/ŝ scaling, the signal is almost two orders above the SM background, as seen in the left panel of Fig. 1. For comparison, in the case of the cross section for the γγ channel shown in the right panel, the difference between signal and background is close to three orders of magnitude. This indicates that the gg → Zγ and γγ signals should be readily detectable at high-energy hadron colliders if gQGC operators appear at the TeV scale. Details of a more complete study are given below. Fig. 1 also shows the next-to-leading order (NLO) SM background due to gg → Zγ via box diagram. Its cross section has been calculated using MadGraph@NLO [42] to calculate the loop diagram automatically. Due to the loop factor, this gg background is suppressed by roughly four orders of magnitude compared to its qq counterpart. In addition to its own contribution, this loop-induced gg → Zγ scattering amplitude can also interfere with the dimension-8 operator amplitude. However, for this interference term to be comparable with the qq background, the dimension-8 operator amplitude would need to be two orders larger than the qq one. In this case, the signal is much stronger than background. In the case that the signal is comparable with the background, i.e., the dimension-8 operator contribution is comparable with the qq one, the interference term is two orders smaller. In either case, the loop-induced SM gg background is negligible. Hence, in the following discussion we focus on the SM qq background and the gg → Zγ signal.
We note that dimension-6 operators could in principle contribute to the same final state pp → Zγ. It could well be that such operators are absent, but their presence would not invalidate our analysis. We recall that the amplitude for a 2 → 2 process is dimensionless. With a coefficient that is ∝ 1/Λ 2 , where Λ is the corresponding new-physics scale, the contribution to the amplitude should scale as s/Λ 2 due to dimensional analysis. Comparing with the dimension-8 contribution that scales with s 2 /Λ 4 , the dimension-6 operator is suppressed in the high-energy signal region √ŝ 1 TeV. Hence dimension-6 operators are not a problem for our study of the collider sensitivities to dimension-8 gQGC operators in this paper. We note that this argument based on dimensional analysis applies to any dimension-6 operator, and is independent of its concrete form.

Angular Distribution and Background Suppression
The SM background mainly comes from t-channel qq → Zγ scattering. This process is intrinsically suppressed at large angles by the intermediator quark propagator, which is ∝ 1/t in the massless limit. For a fixed ratio between the Mandelstam variablesŝ andt, the leading-order cross section for the SM background is given by where δ ∼ 0 is an infrared regulator that scales as 1/ŝ in the high-energy limit. The ratio between the quark mass m q and m Zγ , m q /m Zγ , would in principle provide a natural cut-off value for the scattering angle, but in practice m q /m Zγ is much smaller than typical kinematic cuts for event selection, and hence is irrelevant for our purposes. Although the SM background decreases as 1/ŝ, it will still dominate the event rate if the cut-off scales M i of the gQGC operators are much higher than the 1 TeV used in Fig. 1. For this reason, it is desirable to find kinematic cuts to suppress the SM background.
One such cut-off is provided by the scattering angle in the center-of-mass frame [7], since the Z and photon are emitted as initial state radiation in the t-channel exchange process and hence are typically quite forward. The differential cross section dσ/d cos ϑ [27] is where T q 3 and Q q e are the quark isospin and electric charges. As shown by the black curve in Fig. 2, the SM background indeed peaks in the forward and backward directions, where cos ϑ = ±1. The comparison between the left panel ( √ŝ = 100 GeV) and the right panel ( √ŝ = 1 TeV) shows that higher energy leads to more forward scattering. This property applies to both the Zγ and γγ channels.
For comparison, the signal distribution is much less forward-backward peaked than the SM background. The differential cross sections for the gg → Zγ process generated by the dimension-8 operators are (3.5) These angular distributions are the same as for the gg → γγ case [7] in the limit x Z → 0, reflecting the fact the limit x Z → 0 corresponds to a massless Z boson that is no different from a massless photon.  Since the Z boson mass is significantly smaller than the typical invariant mass m Zγ at hadron colliders, the correction due to a finite x Z is generally negligible.
The differences between the gQGC signals (colorful non-solid lines) and the SM background (black solid line) in Fig. 2 are prominent. Since the polarizations of the incoming gluons are unknown, the process is symmetric under the interchange g( p 1 ) ↔ g( p 2 ). This means that there is a symmetry in the polar angle under ϑ → π − ϑ, under which sin ϑ → sin ϑ and cos ϑ → − cos ϑ. Fig. 2 shows that O gT,(0,4) yield isotropic distributions, whereas the other gQGC operators yield distributions that have small forward-backward peaks. The four Lorentz structures have quite different distributions in general, but we note the following similarities. At the low energy √ŝ = 100 GeV in the left panel, O gT,(0,4) and O gT, (3,7) yield similar angular distributions, whereas in the right panel at √ŝ = 1 TeV there is greater similarity between the angular distributions for O gT, (2,6) and O gT, (3,7) .
Since the gQGC signals have only mild anisotropies, whereas the SM background is concentrated in the forward and backward regions, a simple cut on the Z/γ polar scattering angle in the rest frame of the gg system can suppress significantly the SM background and thereby enhance the signal sensitivity. This feature is used in our study of the experimental sensitivities at the LHC and future colliders in Sec. 4. The potential for distinguishing different gQGC operators via the angular distributions is discussed further in Sec. 6.

Search Strategy and Sensitivities at Hadron Colliders
As discussed above, probes of the gQGC operators at hadron colliders with the processes gg → Zγ and γγ are potentially interesting. The signals would dominate over the SM background when the gg center-of-mass energy reaches the TeV scale, and a simple cut on the scattering angle can further enhance the sensitivity significantly. In this section, we study more details of possible experimental probes. The photon signal at a hadron collider is very clear, so we simply use photon information to obtain sensitivities for the gg → γγ channel. The event spectrum for the gg → γγ process has been studied carefully in [7], and we discuss the updated sensitivity of gg → γγ in Sec. 4.1.  (1,5) O gT, (2,6) O gT,(3,7) The normalized event spectra for gg → Zγ due to the various dimension-8 operators at the LHC with √ s = 13 TeV are shown in Fig. 3. The major differences between Fig. 1 and Fig. 3 are due to the gluon parton distribution function (PDF) in proton, which suppresses the high-energy tails of the event spectra. This makes the difference between the gQGC signals and the SM background even more prominent, and reduces the sensitivity to the treatment of the unitarity constraint. Spectra in the Zγ invariant mass m Zγ are shown in the left panel of Fig. 3 and the photon transverse momentum p T γ distributions are shown in the right panel. 1 Comparing the two panels of Fig. 3, we see that the p T γ spectra grow more slowly. In the cases of charged-fermion final states, Z → + − and qq, the invariant mass spectrum contains more information, by virtue of the p 8 ∼ s 4 energy dependence. However, in the case of the decays into neutrinos, Z → νν, the only reliable experimental information is provided by the photon transverse momentum p T γ .
We study individually the Z boson decays into various fermion-antifermion pairs: decays into charged leptons, Z → + − , in Sec. 4.2; decays into neutrinos, Z →νν, in Sec. 4.3; and decays into quarks, Z → qq, in Sec. 4.4. The prospects for all these channels at the future hadron colliders (HE-LHC, FCC-hh, and SppC) are studied in Sec. 5. The event rates in this paper are simulated with MadGraph5 [28] with model files prepared by Feynrules [29] in the UFO format [30] for the leading order (LO) results. NLO corrections to the background at √ s = 13 TeV are taken into account by an overall K-factor that is extracted by comparing our LO calculations with the NLO calculations given in the experimental papers (the pp → γγ process reported in [33], the pp → Zγ processes in [20] for the leptonic channel and [21] for the invisible decay channel as well as [22] for the hadronic channel). For future colliders, we have used the same overall K-factor.
We adopt the same binning for the prospective signals, N sig i , as has been used for the corresponding experimental data points, N data i . On the basis of these numbers, we use the following χ 2 function to evaluate the sensitivity to the dimension-8 gQGC operators: 1 Analogous mγγ distributions for gg → γγ scattering were shown in Fig. 2 of [7].
In addition to the cut-off scale M i that enters implicitly through the signal event rates N sig (M i ), the χ 2 function also contains an overall normalization F and multiple nuisance parameters f k to account for the various systematics published in the experimental papers. For each systematical error, both the central valueσ i and the deviations f k σ k,i have been taken into account. Assuming quadratic dependences on the nuisance parameters F and f k , the χ 2 minimization can be done analytically [31,32]. After marginalizing over F and f k , the resulting χ 2 (M i ) as function of a single parameter provides the sensitivity to the cut-off scale M i . The same procedure has been followed for the SM background, so as to calculate χ 2 min . The square root of the difference, ∆χ 2 (M i ) ≡ χ 2 (M i ) − χ 2 min , directly gives the significance that we show below. We generally quote results at the 95% C.L., corresponding to ∆χ 2 = 3.84. The sensitivities shown in this paper are obtained for each of the 8 gQGC operators separately and for the BISM combination. The constraints on the individual operator coefficients would in general be weaker if all the operators were included simultaneously and one marginalized over the other operator coefficients.
Below we apply this χ 2 analysis to the ATLAS data for various channels. To guarantee that the quadratic χ 2 function is a good enough approximation, only those bins with at least 5 events are included in Eq. (4.1). This provides a very conservative sensitivity estimation, since beyond the end points of the SM background (or the experimental data), the signal is still growing and essentially free of background. If a more sophisticated data analysis were employed, the sensitivity could be further improved. We leave this to our experimental colleagues.

pp → γγ at the LHC
We first update the analysis of [7] using the ATLAS measurement of the diphoton invariant mass m γγ spectrum using 36.7 fb −1 to include 139 fb −1 [33], which provides improved sensitivities to dimension-  [34]. The constraints at the same 95% C.L. can improve by amounts ∼ 1 to 2 TeV.
The solid black line in Fig. 4 represents the sensitivity to the Born-Infeld O gT,BI combination of the gQGC operators in Eq. (2.5). This linear combination of 6 individual gQGC operators yields an event rate that is much larger than any individual operator, as seen in Fig. 1. Consequently, the 95% CL lower limit on the BISM mass scale reaches M BI 5 TeV. This result is strong enough to impose a significant constraint on the separation between branes in some M-theory inspired models that address the electroweak hierarchy problem.
We emphasize that we are working to quadratic order in the dimension-8 operator coefficients, unlike many analyses of dimension-6 operators at the LHC. Also, our analysis is reliant on values of √ŝ that are below those where the unitarity constraints become important. However, we also note that in some instances the values of √ŝ contributing to our analysis are comparable to (or even exceed) the magnitudes of the constraints on new-physics scale M i that we derive. In such cases, the validity of   The CMS and TOTEM collaborations have recently searched for exclusive diphoton production in proton-proton collisions and set constraints on the dimension-8 quartic photon couplings ζ 1 F µν F µν F ρσ F ρσ and ζ 2 F µν F µρ F ρσ F σν [18]. Their result is a two-dimensional contour in the plane of the coefficients ζ 1 and ζ 2 , in the absence of theoretical constraints. As demonstrated earlier, the Born-Infeld QED extension in Eq. (1.1) naturally correlates the two purely photon QGC operators with ζ 1 = ζ 2 = 1/32M 4 . Thus the 95% sensitivity of the CMS-TOTEM data becomes M BI ≥ 670 GeV. The CMS and TOTEM sensitivity is limited by the fact that it uses γγ → γγ scattering with photons in both the initial and final states, leading to a weaker constraint than the gg → γγ process at the same collider for either individual operator or the combined BISM operator. This confirms the advantage of considering the gg → γγ channel at hadron colliders.

pp → Z( + − )γ at the LHC
The ATLAS Collaboration has searched for new physics in the Z( + − )γ final state with an integrated luminosity of 139 fb −1 [20] at √ s = 13 TeV. Here we re-interpret the results with electron (Z → e + e − ) and muon (Z → µ + µ − ) final states as constraints on the gQGC operators. Both the photon and the charged leptons are observable. The event selection required the photon to have pseudo-rapidity in the range |η γ | < 2.37 and transverse energies E γ T > 30 GeV. We use the ATLAS data with tight photon identification, which had an identification efficiency ranging from 82 − 85% for E γ T ≈ 30 GeV to 90 − 98% for E γ T > 100 GeV. The electrons and muons are required to have pseudo-rapidities |η | < 2.47 and transverse momenta p T > 25 GeV. The identification efficiency for charged leptons is about 80% (93%) for p T ≈ 25 GeV (100 GeV). Because of geometrical limitations, the transition region between the barrel and endcap (1.37 < |η| < 1.52) is excluded for both the photon and the charged leptons. The dominant instrumental backgrounds for the Z( + − )γ final states are Z + jets, pile-up events with one photon, ttγ and τ + τ − γ, as well as double vector boson production with or without an isolated photon [20]. It was estimated that the combination of these backgrounds contributes about 18% of the signal process qq → Z( + − )γ in the fiducial phase space. For simplicity, we take the Sherpa simulation of the background from the ATLAS paper [20] and assign an overall normalization factor F when fitting the experimental data as explained in the discussion of Eq. (4.1).
The left panel of Fig. 5 shows the sensitivities obtained from fitting the observed Z(→ + − )γ invariant mass m Zγ spectrum at ATLAS [20] for probing the gQGC operators. Since the contributions of GG i W i W i and GGBB operators with the same Lorentz structures have identical EW mixing factors for the Zγ combination, except for a relative sign, see Eq. With three detectable particles in the final state, much more information can in principle be extracted beyond just the invariant mass, pseudo-rapidity, and transverse momentum. For example, ATLAS also analyzed the angular distribution of the + − pair relative to the scattering plane. We discuss in Sec. 6 potential improvements of the cross-section analysis using spin and angular correlation measurements.

pp → Z(νν)γ at the LHC
The detection of the Z(νν)γ final state differs from those of the diphoton channel and the gg → Z( + − )γ mode, since only the photon is observable. This renders impossible the reconstruction of the invariant mass spectrum that is the optimal choice for seeing the momentum dependence of 2 We recall that the Zγ final state does not provide any constraint on the Born-Infeld scale, since the contributions of the W 3 W 3 and BB operators cancel with each other as explained below Eq. (2.2). the dimension-8 gQGC operators. The only usable information is the photon energy/momentum vector. Although this can readily be measured, the information cannot be fully used since the initial momenta of the colliding gluons in the beam protons are unknown. Because of this uncertainty, one can only utilize the transverse momentum spectrum depicted in the right panel of Fig. 3. Although not optimal, it still possesses the feature that the event rate keeps growing with momentum until the unitarity saturation is encountered and the gluon PDF suppression finally dominates.
The ATLAS measurement of Z(→νν) + γ uses an integrated luminosity of L = 36.1 fb −1 at 13 TeV [21]. The event selection uses the following experimental cut on the transverse energy: E γ T , E miss T > 150 GeV. Similarly to the charged lepton case, the photon pseudo-rapidity is in the range |η γ | < 2.37, but excluding 1.37 < |η γ | < 1.52 to avoid the gap between the detector barrel and end-cap regions. The transverse missing momentum p miss T needs not be opposite to the photon transverse momentum p γ T , due to parton shower from the initial gluons. Requiring the azimuthal angular difference ∆φ(p γ , p miss T ) > π/2 further suppresses the background.
The dominant instrumental background for the Z(νν)γ final state comes from W γ associated production with the W decaying leptonically but the charged lepton from W → ν not being identified. In addition, jets +γ and other events with either an electron or a jet being misidentified as a γ are possible sources of background [21]. These backgrounds can amount to about 68 (58)% of the inclusive (exclusive) SM background qq → Z(νν)γ in the fiducial phase space, i.e., the cross section with ≥ 0 (0) additional jets.
The sensitivities from analytical χ 2 fits to the different gQGCs obtained using Eq. (4.1) are shown in the left panel of Fig. 6. Although the pp → Z(νν)γ channel is not the optimal one, the lower bound on the cut-off scale M at 95% C.L. can still reach around 1.4 TeV for O gT,(1,2,3,5,6,7) and 1.9 TeV for O gT,(0,4) . One reason is that the Z →νν mode includes all three neutrinos whereas only electron and muon modes are used for the charged lepton mode. Moreover, the Z decay branching ratio into a single neutrino is roughly double that of its charged lepton counterpart, so the combined Z →νν branching ratio is about three times larger than the combined Z → + − branching ratio. The right panel of Fig. 6 shows the prospects at the HL-LHC with its higher energy and luminosity. The lower bounds are improved by 1 TeV.

pp → Z(qq)γ at the LHC
The leptonic modes discussed in the previous sections have very clean signals at the hadron colliders. However, although the hadronic decay modes of the Z boson have more complicated backgrounds, the larger branching ratio provides some advantages. The ATLAS Collaboration has measured the Z(→qq)γ m Zγ invariant mass spectrum with an integrated luminosity of L = 36.1 fb −1 at 13 TeV [22]. The experimental analysis first selected events containing a hadronic jet with p T > 200 GeV and |η| < 2.0 in addition to a photon with p γ T > 250 GeV and |η| < 1.37. The events were then divided into four categories [22]: 1) The BTAG category of events in which the two leading track-jets associated with a large-radius jet candidate satisfy the jet mass requirement for b-tagging, with a Z →qq identification efficiency of 3 to 4%; 2) The D2 category that is composed of events satisfying combined jet mass and substructure discriminant requirements with an identification efficiency of 20 to 28%; 3) The VMass category containing events that pass the jet mass selection but fail to enter either of the two previous categories (identification efficiency 24 to 36%); and 4) The Else category containing the remaining events that pass the baseline selection (identification efficiency 40 to 50%). Roughly speaking, the event rates in the above four categories increase by an order of magnitude going from category 1) to category 4). Fig. 7 shows some validation results from our simulation obtained using Pythia8 [35] for parton shower effects and Delphes3 [36] for detector effects, as well as FastJet [37], ExRootAnalysis [28] and Root [38] for data analysis. The invariant mass distribution in the left panel reproduces adequately the Z boson peak and a lower background peak around 20 GeV. Since the event rate is dominated by the VMass and Else categories, we only show these two in the middle and right panels. Our simulations are reasonably consistent with the ATLAS data points.
The sensitivities for the pp → Z(qq)γ channel at the LHC are shown in Fig. 8 for both the current LHC data (left panel) and the future HL-LHC (right panel). Although the hadronic modes are much more difficult to detect, the 95% C.L. lower bound still reaches around 1.5 TeV for O gT,(1,2,3,5,6,7) and 2.3 TeV for O gT,(0,4) , which are even stronger than the neutrino mode. Further improvement by at least another 1 TeV should be possible at HL-LHC. In addition to the pseudo-rapidity and transverse momentum, the angular distribution of theqq pair relative to the scattering plane can also be measured. However, this was not done in the ATLAS analysis [22]. A discussion of possible future improvements using this information is given in Sec. 6.

Prospective Sensitivities at Future Hadron Colliders
In addition to the existing LHC and its upgrade to HL-LHC, other hadronic colliders are being proposed, including HE-LHC [34], FCC-hh [39], and SppC [40], whose collision energies range from 25 TeV to 50 TeV and 100 TeV. As discussed in Sec. 3.1, because of its strong momentum dependence, the gQGC signal increases very rapidly with the colliding energy. This gives future higher-energy hadron colliders great advantages for probing the gQGC operators. In addition, future hadron colliders typically have much larger luminosity than the current LHC. We assume a universal figure of L = 20 ab −1 for the various higher-energy colliders and higher p T cutoffs than for the LHC, as shown in Table 1, but retain the same scattering-angle cuts as for the LHC.     Fig. 9: The 95% C.L. sensitivities of the various channels at the LHC, HL-LHC and proposed future hadron colliders, with cuts given in Table 1.
the different colliders are colour-coded with shadings to distinguish the various final states considered. As expected, the higher the collision energy, the greater the sensitivity. The greatest Zγ sensitivity is in general provided by Z → νν decays, while the γγ channel provides even greater sensitivities for O gT, (4,5,6,7) .
We also show the combined sensitivities of the Zγ and γγ channels by summing the χ 2 functions in Fig. 4, Fig. 5, Fig. 6, and Fig. 8 for the LHC and similarly for the future higher-energy colliders. The lower bound at 95% C.L. can be significantly enhanced in 100 TeV collisions to above 20 TeV for O gT,(0,4) and the SM Born-Infeld extension, an order of magnitude beyond the current LHC limits.

Potential Improvements using the Z Spin and Angular Correlation
As mentioned above, the pp → Z( + − ,qq)γ channels have three detectable particles in the final state. In these cases, not only the scattering angle distribution at the intermediate Zγ level that is shown in Fig. 2 can be used for background discrimination, but also the difference between the Zγ scattering plane and the decay plane of Z → + − , qq can provide useful information. Thus the phenomenology of the Z( + − ,qq) channels is in principle much richer than that for pp → γγ and pp → Z(νν)γ, which yield only photons in the final state. In this section we explore basic features of the Z spin and angular correlations with a view of potential improvements in the analysis sensitivity.
In order to keep full information on the Z → + − , qq decays, we adopt the helicity amplitude formalism [41] to calculate the scattering matrix element for the gg → Zγ process: Here, λ 1,2 = ±1 and λ γ = ±1 are the helicities of the massless gluons and photon, respectively, whereas in the case of the massive Z boson we have λ Z = 0, ±1. Also, λ f /2 = ±1/2 denote the helicities of the final-state fermions from Z decay. For simplicity, we assume that the masses of the final-state fermions can be neglected.
The amplitude for the two-step process of Zγ production followed by Z boson decay can be written as where P ν ≡ Z * µ γ|O i |g 1 g 2 and D µ ≡ ff |O SM |Z * ν are the matrix elements for the production and decay of Z * , respectively. The intermediate Z-boson propagator G µν (Z * ) can be written as Since the decay width of the Z boson is considerably smaller than its mass, we use here the narrowwidth approximation ŝ − m 2 Z + imΓ Z G µν (Z * ) ≈ λZ µ λZ (p) ν * λZ (p). Then one polarization vector ν * λZ contracts with P ν to give the production matrix element while the other polarization vector µ λZ contracts with D µ to give the decay matrix element. Consequently, the total helicity amplitude (6.2) can be decomposed into a sum of products of three scalars: where D Z ≡ŝ − m 2 Z + im Z Γ Z . The summation is over the helicity of the internal Z boson, which does not change during free propagation. The production and decay matrix elements are which can be calculated in the usual way for gg → Zγ and Z → ff separately.
The spin information appears in the absolute square of the scattering matrix element, The Z-boson polarization vectors encode the spin information.
Resolving the decay plane of the channel Z →νν is not possible, and resolving that of Z →qq is challenging in practice, particularly when the Z boson is highly boosted. Therefore here we investigate Fig. 10: Kinematic variables for the process gg → Z(ff )γ. The scattering plane is defined by the incoming gluons and outgoing Z/γ bosons. The polar angle (scattering angle) of the Z boson is defined as ϑ, z * is defined as theẑ * -axis in the rest frame of the Z boson along the direction of the Z boson momentum, and thex * −ẑ * plane is spanned by the vector z * and the direction of the initial gluon motion. The polar and azimuthal angles of the leptons in the rest frame of the Z boson are θ * and φ * .
in detail only the leptonic decay channel. The Z-boson momentum provides a natural definition of the z-axis with respect to which we define the fermion polar angle θ in the Z rest frame. Since the detector is azimuthally symmetric, we can without loss of generality define thex-axis to lie in the scattering plane of the Z and γ bosons, with the Z momentum having a positivex component. The azimuthal angle difference between the scattering and decay planes is φ , see Fig. 10. Then the decay polarization tensor components D λ Z λZ are The polarization state of the intermediate Z boson depends in general on the dimension-8 gQGC operators that produces it, carrying information beyond that provided by its momentum. 3 The decays of the Z boson open a window on the polarization information. Below we show how to use the fermion polar angle in Sec. 6.1 and the azimuthal angle in Sec. 6.2 to distinguish the gQGC operators. With enough event rate, it is also possible to study a combined two-dimensional distribution of polar and azimuthal angles as we discuss in Sec. 6.3.

Polarization Effects in the Fermion Polar Angle Distribution
The Z polarization tensor P λ Z λZ is defined in the laboratory frame with the z-axis along the beam direction. Since the scattering is rotationally-invariant with respect to the azimuthal angle, the components of P

256
(1 − x Z ) 2 16 x Z sin 2 ϑ 5 + 3 cos(2ϑ) , (6.9h) for the 8 gQGC operators. We see that the positive (λ Z = +) and negative (λ Z = −) polarization states are always produced equally. In other words, there is no net polarization effect. Furthermore, the production rate of the longitudinal state (λ Z = 0) is suppressed by x Z at large invariant mass m Zγ , and therefore can be neglected at high-energy hadron colliders. However, in the region x Z ∼ O(1) it could be useful for distinguishing different Lorentz structures.
In the presence of the longitudinal mode, there is in principle a polarization effect, though its magnitude decreases with increasing center-of-mass energy m Zγ . The polarization of the Z boson can be measured by studying angular distributions of its decay products. The differential cross sections with respect to the polar angle of the outgoing lepton are given by dσ γ + − ,(0,4) σ γ + − ,0,4 d cos θ = 3 8 1 + cos 2 θ , (6.10a) dσ γ + − ,(1,5) σ γ + − ,1,5 d cos θ = 3 26(1 + cos 2 θ ) + 16x Z sin 2 θ + 3x 2 Z cos(2θ ) + 9x 2 For O gT,(0,4) , since the Z boson is generated isotropically as shown in Fig. 2, the polar angle of the lepton is identical with the distribution in the decay of an unpolarized Z boson. However, the situations for other operators differ because of the contributions of the longitudinal component. The left panel of Fig. 11 shows the normalized differential cross sections at the gluon-gluon collision energy √ŝ = 100 GeV, where sizeable differences can appear. However, these differences are proportional to the scale factor x Z , and in the limit x Z → 0 all the curves in the left panel of Fig. 11 converge to forms ∝ 1 + cos 2 θ .

Spin Correlation in Fermion Azimuthal Angle Distributions
Since the intermediate Z is not directly measured, it is in general produced in an entangled spin state with the only exception being when produced by the operators O gT,(0,4) . Hence different polarizations correlates with each other and the interference between the helicity states can induce non-trivial distributions of the fermion azimuthal angle. To make this transparent, it is necessary to calculate the relevant off-diagonal matrix elements of the gg → Zγ production process: The related off-diagonal matrix elements of the Z → + − decay process are given in (6.8c) and (6.8d). We see that the transverse-transverse correlation effect leads to a non-trivial dependence on the azimuthal angle φ of the outgoing lepton. The normalized differential cross sections with respect to the azimuthal angle of the outgoing lepton for the different operators are A non-trivial distribution of the azimuthal angle φ * would indicate the importance of a gQGC operator different from O gT,(0,4) , offering in principle the possibility of distinguishing the gQGC operators. The right panel of Fig. 11 shows the normalized differential cross sections at a gluon-gluon collision energy √ŝ = 100 GeV. We see that, compared to the polarization effect in the left panel of Fig. 11, the correlation effect is expected to be more sensitive to the Lorentz structures of the gQGC operators. However, the effect is suppressed by x Z at high M Zγ .

Conclusions
The SMEFT is a powerful way of searching for possible indirect effects of new physics beyond the SM that may appear at energies outside the direct reaches of active experiments. In particular, the SMEFT provides a framework for exploring possible anomalous multi-boson couplings while taking into account all the gauge symmetries of the SM, which is one of the priority measurements for LHC and possible future colliders. The present and prospective experimental sensitivities to anomalous triple gauge couplings (TGCs) and quartic gauge couplings (QGCs) have been studied extensively, but there have been fewer studies of possible quartic couplings between gluons and electroweak gauge bosons, which are absent in the SM. In the SMEFT, these are first generated by dimension-8 operators, whose study offers an interesting window on BSM physics [7], complementing the exploration of dimension-6 operators that have been analyzed extensively, see [4,5] and references therein.
In this paper we have presented a first analysis of the experimental sensitivities of measurements of the gg → Zγ process to the dimension-8 quartic couplings of pairs of gluons to the Z and photon (gQGCs), including an analysis of present LHC data and assessments of the prospects at HL-LHC and proposed higher-energy proton-proton colliders. Four distinct Lorentz structures of CP-conserving dimension-8 operators contain gQGCs, and each may be constructed either with pairs of SU(2) or U(1) gauge field strengths. We have stressed the differences between the polar angle distributions they generate from the one generated by the dominant SM background due toqq → Zγ, and analyzed the possible sensitivities of analyses of Z → + − ,νν andqq final states. The present LHC data on these final states correspond to integrated luminosities up to 139 fb −1 at 13 TeV, and we show that they exclude mass scales 2 TeV in the dimension-8 operator coefficients, with future colliders offering sensitivities to scales an order of magnitude higher, extending above 20 TeV for O gT,(0,4) . The gg → ZZ channel has more complex phenomenology that will be studied in a separate paper.
We also present an updated analysis of the sensitivities of LHC measurements of the γγ final state with up to 139 fb −1 at 13 TeV and at future colliders. The present data already constrain the nonlinearity scale of the Born-Infeld extension of the SM to be 5 TeV, putting pressure on some brane models that address the electroweak hierarchy problem. A future collider with 100 TeV in the center of mass could be sensitive to Born-Infeld scales 20 TeV. As we also show, potential analysis improvements might be made possible by exploiting Z spin measurements and angular correlations using decays into + − andqq final states, which are important at low Zγ invariant masses.
Our results provide further illustrations of possibilities for exploring and constraining possible dimension-8 terms in the SMEFT, in addition to light-by-light scattering [6,18] and the searches for neutral triple-gauge couplings discussed elsewhere [8,9]. There is plenty of life in the SMEFT beyond dimension 6.