Quartic Gauge-Higgs couplings: Constraints and Future Directions

Constraints on quartic interactions of the Higgs boson with gauge bosons have been obtained by the experimental LHC collaborations focussing on the so-called $\kappa$ framework of flat rescalings of SM-like interactions in weak boson fusion (WBF) Higgs pair production. While such approaches are admissible to obtain a qualitative picture of consistency with the SM when the statistical yield is low, once more statistics become available a more theoretically consistent framework of limit setting is desirable. Reviewing the constraints provided at the Large Hadron Collider, we first show that these limits are robust when considered in a leading order context. Turning to radiative corrections, we demonstrate the limitations of this approach in the SM, and by adopting Higgs effective field theory techniques, we clarify the sensitivity from single Higgs measurements to rescalings of quartic Higgs-gauge couplings. We then discuss avenues for sensitivity improvements of WBF analyses employing Graph Neural Networks to combat the large contributing backgrounds.


Introduction
Ten years into the Higgs characterisation programme at the Large Hadron Collider (LHC), many couplings of the Higgs boson to other SM matter have been shown to largely follow the Standard Model (SM) expectation. Given the insights from electroweak fits before the Higgs discovery in 2012, this was perhaps anticipated for the gauge sector [1,2]. The observation of H → γγ with a rate in agreement with the SM prediction indicated alignment of the fermion-Higgs interactions as part of the Higgs discovery for the top quark. Later analyses in other fermion decay modes of the Higgs boson [3][4][5][6] have further established the Higgs boson as SM-like.
Couplings that have received less attention, also because corresponding searches are statistics-limited, are the quartic interactions of the Higgs boson with gauge fields. Assuming a weak doublet-like character of the Higgs boson, these interactions are fully correlated with trilinear Higgs-gauge (HW W , HZZ) interactions due to gauge invariance in the SM. Searches for anomalous quartic interactions in the so-called κ framework [7] are therefore theoretically cumbersome, albeit instructive when the statistical sensitivity is low (similar to κ-like analyses of single Higgs observables after the Higgs boson's discovery).
In the expansion of the Higgs boson around its vacuum expectation value (vev), the gauge boson masses are a property of the non-linear SU (2) L × U (1) Y electroweak symmetry realisation, whilst the Higgs bosons' couplings to W and Z bosons probe the alignment of quantum fluctuations around this vev. These interactions carry important information as any departure from the SM expectation immediately implies perturbative unitarity violation at a scale above the vev [8][9][10], which acts as a strong indicator of a new scale of physics beyond the SM (BSM).
Quartic interactions ∼ H 2 V 2 (V = W, Z), while being less relevant from a perturbative unitarity perspective, still carry important information when HV V couplings (dis)agree with the SM: Different BSM scenarios show similar alignment of single Higgs interactions with the vev. In the gauge sector this degeneracy is only lifted by considering quartic interactions. An example for this is custodial singlet mixing [11][12][13][14] compared to minimal compositeness Higgs models (MCHMs). In the former case Higgs interactions are modified by a characteristic mixing angle ∼ cos θ whereas in, e.g., MCHM5 [15], we obtain a modification of ∼ 1 − v 2 /f 2 [15][16][17] (where v 246 GeV is the vev and f v is the decay constant of the Callen-Coleman-Wess-Zumino (CCWZ) [18,19] construction). On the sole basis of a precision analysis of the HW W, HZZ interactions, it is therefore impossible to discriminate between Higgs mixing and Higgs compositeness as we can always understand cos 2 θ = 1 − v 2 /f 2 (note that custodial symmetry is sufficiently conserved in both models [20]). Only when we consider the curvature of the Higgs interactions is this ambiguity lifted in the gauge sector: in the custodial singlet mixing scenario, the H 2 V 2 vertex is modified by cos 2 θ = 1 − 2v 2 /f 2 , where the latter is the result in MCHM5.
Targeted analyses for flat rescalings in weak boson fusion (WBF) Higgs pair production show that this channel will be statisticslimited at the high-luminosity (HL) frontier [21][22][23][24][25]. To some degree, this is similar to investigations of the Higgs boson's self-coupling; statistics is also a limiting factor at the HL-LHC for rescalings of the trilinear coupling κ 3 [26]. To partially address this issue, single Higgs probes have been suggested [27][28][29][30][31][32] to add sensitivity. The comparison of single-Higgs probes (or even electroweak precision data [33,34]) to rescalings of the Higgs self-coupling κ 3 should, however, also be understood as a theoretical figure of merit related to the trustworthiness of the obtained constraints as different orders in the perturbative (effective field theory) series expansion are tensioned against each other. Given that WBF Higgs pair production is significantly smaller than inclusive Higgs pair production via gluon fusion [35], similar issues of κ 2V do occur. It is the purpose of this work to clarify these questions in the light of existing analyses pursued by the ATLAS and CMS experiments [36][37][38]. We show that this is achieved by understanding the SM as a particular parameter point of the electroweak chiral (or Higgs Effective Theory, HEFT) Lagrangian. This enables us to consistently modify the quartic Higgs-gauge interactions away from the SM expectation with theoretically consistent implications for single Higgs observables. In particular, we will focus on radiative corrections to the H → V V rates, which we tension against direct sensitivity obtained from WBF Higgs pair analyses. To formulate an as-good-as-possible direct search limit, we employ machine learning techniques to the signal-background discrimination to show that significant improvements are possible beyond the existing, more traditional techniques discussed in the literature. This work is organised as follows: In Sec. 2, we review the theoretical baseline of κ 2V analyses in the SM. We consider electroweak precision constraints from oblique corrections and derive unitarity constraints from a partial wave analysis, clarifying the validity of the results of [36][37][38] viewed as a leading order signal analysis. We show that a κ 2V analysis cannot be theoretically defended when we consider quantum corrections in the SM. To overcome this limitation, we turn to the electroweak chiral Lagrangian in Sec. 3, which provides the theoretically consistent playground to discuss correlations between single and double Higgs modifiers at a given order of the amplitudes' loop expansion. Performing the next-to-leading order calculation of H → ZZ * , W W * alongside H → γγ, γZ, we discuss the current and expected sensitivity of single Higgs data to κ 2V coupling modifiers. This highlights the direct sensitivity provided in weak boson fusion HHjj analysis as the sensitive tool to constrain κ 2V in the future (for a recent discussion at e + e − machines see Ref. [39]). Any improvement beyond increasing statistics of current strategies [36][37][38] will therefore be directly correlated with an improved understanding of the Higgs boson's TeV scale characteristics. To this end, we show that the application of state-of-the-art Graph Neural Network (GNN) techniques [40,41] to pp → HHjj final states can greatly reduce the dominant backgrounds of such searches, thus leading to a significant sensitivity enhancement to κ 2V . We summarise and conclude in Sec. 6.

Vices and Virtues of κ 2V
We consider first the effect of flatly rescaled quartic Higgs-gauge couplings to electroweak precision constraints. Here, we particularly focus on oblique corrections [42,43] and their sensitivity to coupling modifications of only the physical Higgs boson (with mass M H 125 GeV). These interactions are not associated with momentum dependencies at one-loop order so we directly obtain ∆S = ∆U = 0 , (2.1) (we do not consider single Higgs modifiers at this point). The T parameter is sensitive to the custodial isospin violation for κ 2W = κ 2Z with Λ as the UV cutoff regulator for demonstration purposes. 1 This means that for large cutoffs of order 10 TeV we need to ensure κ 2W κ 2Z at the 1.5% level in order to not violate present constraints [44] whilst UV engineering interactions that dominantly source custodial isospin violation through the quartic gauge-Higgs interactions. The latter is naively possible in the κ framework, Eq. (2.2), but rather unlikely when considered in the context of more realistic SM extensions, some of which we have alluded to above. On the other hand, by identifying κ 2W = κ 2Z = κ 2V , constraints from oblique corrections are removed at the  considered one-loop level. This is the identification we will use in the following; electroweak precision results suggest that the Higgs boson should be considered as a custodial iso-singlet state, which is the hypothesis on which the electroweak chiral Lagrangian is built. Next, we turn to unitarity constraints. The partial waves of i 1 i 2 → f 1 f 2 scattering process for angular momentum J are given by (see Ref. [45]) (modulo factors of 1/ √ 2 for identical particles in the initial or final state) where we denote √ s as the centre-of-mass energy, the scattering angle in the centre-of-mass system is θ, and The D J µ i µ f are the Wigner functions of the Jacob-Wick expansion [45] related to the helicities of the initial and final states Unitarity of the S matrix then translates for f = i to the familiar conditions The process relevant for κ 2V constraints is longitudinal HV L → HV L scattering, so that the Wigner functions reduce to the Legendre polynomials, with the J = 0 partial wave providing the dominant constraint. This is shown in Fig. 1, where we define We also include the LHC results from the ATLAS and CMS experiments of Refs. [36][37][38].
As can be seen, the unitarity constraints are relatively loose. The measurements, taken as a leading order result, can be considered perturbatively consistent as the energy scales explored by these analyses do not probe the unitarity cutoff of the leading order theory. The discussion so far suggests that κ 2V is indeed relatively unconstrained when naively considered as a free parameter in SM-like analyses. However, beyond tree-level, a choice of κ 2V = 1 leads to incurable shortfalls within the SM. For instance, when considering electroweak corrections to H → ZZ * in general R ξ gauge [46], we obtain a counter term amplitude (M CT )-renormalised 2 virtual amplitude (neglecting fermion contributions and considering M H > 2M Z for convenience) (Z i ) are the transverse polarisations of the Z bosons (related to the decay current Z → ff for massless fermions) and is the usual MS-related divergence in dimensional regularisation D = 4 − 2ε [47]. This result is qualitatively not surprising: As we have modified without much care only a single part of the gauge-invariant SU (2) L × U (1) Y kinetic Higgs field term 2 being the Higgs doublet) rather than the whole term, we have lost gauge invariance entirely. This is highlighted by the non-vanishing ξ Z gauge fixing parameter dependence of the H → ZZ decay. ζ 2 = 0 has also induced singularities that are not mended by a renormalisation of the bare HZZ SM interactions (included in Eq. (2.7)). 3 Again this is expected; radiative ζ 2 corrections should renormalise coupling deviations (the Lorentz structure of Eq. (2.7) is linked to the HZZ interactions of Eq. (2.8)). However, even if we introduce a bare ζ 1 , its renormalisation will be gauge-dependent. Therefore, ζ 1 cannot directly be associated with a physical observable (in contrast to, e.g., the electromagnetic coupling e and its behaviour under renormalisation group flow). The combination of naively weak κ 2V tree-level unitarity constraints and the breaking of gauge invariance within the SM prompts us to a practical phenomenological issue when considering the analyses by ATLAS and CMS (or future experiments) within the κ framework: What is their phenomenological value and what is the level of ζ 2 deviations that can be reasonably expected?
Naturally, this question is not a sensible one in the context of the Standard Model, as highlighted by Eq. (2.7). In concrete UV extensions, some of which we alluded to above, the answer is immediately provided by examining the H 2 V 2 operators and their relation to the theories' input parameters; they will not be free parameters and typically be predictions of BSM electroweak input parameters that are more precisely determined from other measurements. 4 As the analyses and measurements of [36,38] are model-independent in the sense that they modify certain parameters whilst keeping others entirely SM-like, our task is to house these assumptions in a theoretically meaningful way. This requires an approach with which correlations are not plagued by inconsistencies displayed by Eq. (2.7) to inform the reach of κ 2V deviations in the light of constraining single Higgs measurements.

Relevant HEFT interactions at NLO
The only theoretically plausible way to treat H n V V interactions as independent is by avoiding the physical Higgs boson being part of an electroweak doublet, which also means that Standard Model EFT (SMEFT) [48] is not applicable. This is possible by considering the electroweak chiral Lagrangian (or HEFT) [49][50][51][52][53][54], which has received considerable attention recently, including the treatment of radiative corrections [55][56][57][58][59][60][61]. HEFT is intrinsically different from the SM, which becomes clear when radiative corrections are considered in general gauge (chiefly analysed by Herrero and Morales in [60][61][62]), but contains the SM as a particular parameter point choice. This enables us to achieve a theoretically consistent correlation of different Higgs channels including radiative corrections in the modern sense of renormalisability (i.e. renormalisation at a given loop order).
As in any QFT, radiative corrections will source higher-dimensional operators, yet HEFT is non-renormalisable in the classical sense. It is therefore necessary to include all relevant operators in the renormalisation programme from the start as their renormalisation is required for a consistent final one-loop result. The κ V (κ 2V ) analysis that we would like to achieve is then related to a particular choice of couplings for the one-loop renormalised HEFT Lagrangian at a given renormalisation scale, as we will discuss below.

HEFT Lagrangian
Following the notations of Refs. [53,60,61,63], we consider first the leading order (chiral dimension-two [64,65]) HEFT Lagrangian. W a µ and B µ are the SM gauge bosons associated with the SU (2) L × U (1) Y gauge symmetry. Contrary to the SM, the Higgs boson is a singlet field in HEFT and the Goldstone bosons π a are parametrised non-linearly using the matrix In the SM these couplings are rather obviously determined by the choice of precisely measured electroweak input parameters, e.g. gHZZ = eMZ /(sW cW ), gHHZZ = e 2 /(2s 2 where τ a are the Pauli matrices with a = 1, 2, 3. The U matrix transforms under L ∈ SU (2) L , U (1) Y ⊂ SU (2) R R as U → LU R † and is expanded as The gauge fields in the physical basis are defined as The leading order HEFT Lagrangian relevant for this study is given as follows where the interactions of Higgs with gauge and Goldstone bosons are parametrised using the function F H given as ζ 1 and ζ 2 are the new physics parameters (the choice ζ 1,2 = 0 corresponds to the gauge-Higgs interactions of the SM). The ellipses denote the interaction terms with more than two Higgs bosons which are not relevant for our work. L ferm parametrises the fermion-gauge boson interactions, which we take SM-like. L GF and L FP are the gauge-fixing and Faddeev-Popov terms [66], respectively. V (H) is the Higgs potential For our analysis, we will take the potential to be SM-like, κ 3,4 = 1. The gauge-fixing term is written in terms of gauge fixing functions as where the gauge fixing functions are defined as linear with ξ W , ξ Z and ξ A denoting the gauge-fixing parameters in the R ξ gauge. Using these above gauge-fixing functions, the Faddeev-Popov part of the effective, quantised HEFT Lagrangian is then given as where c j are the ghost fields and α j are the corresponding gauge transformation parameters which are similar to that of SM. We parametrise the couplings of Higgs boson with fermions via where y u ij and y d ij are the Yukawa couplings. In this study, we will focus on the case c = 1, which after diagonalisation, directly links the Yukawa couplings to the fermion mass eigenstates with top and bottom masses m t,b , respectively. We will neglect the light quark flavour and lepton masses throughout this work since these do not impact our phenomenological findings.
When considering radiative corrections in the effective theory described by Eq. (3.5), ultraviolet divergences will source new operator structures at one-loop level. This requires the introduction of bare terms to the Lagrangian ab initio for a consistent, one-loop renormalisation programme [61,62]. The terms relevant at NLO one-loop level (chiral dimension 4) are with operators that require UV regularisation listed in Tab. 1.

Amplitudes and Renormalisation
Throughout, we will work in dimensional regularisation D = 4 − 2ε [47] and adopt the on-shell scheme for propagating degrees of freedom [67]. Expanding out the non-linear sigma model of Eq. (3.5), we understand the electroweak vev as a derived quantity with following Sirlin [68,69]. These can be parametrically fixed through choices of e = √ 4πα, M W and M Z . As we identify these quantities as the input parameters for the gauge part of the Lagrangian (coupling constants of SU (2) L × U (1) Y are parametrised as g W = e/s W , g = e/c W , respectively), we have to supply counter terms δZ e , δM 2 bare Lagrangian quantities (we indicate unrenormalised quantities with a subscript '0' in Eq. (3.5), (3.12) after L → L 0 ) e 0 = Z e e = (1 + δZ e )e , The renormalisation of the Weinberg angle is then given by The field renormalisations of the bare Lagrangian quantities are for the photon-Z sector. Tadpoles are renormalised by requiring that they fully cancel against the tadpole counter term. They are treated as part of the parameter renormalisation [67,70,71] so that the tadpole counter term only enters the unphysical Goldstone sector and we can ignore associated contributions in the following. The HEFT coefficients are renormalised through We choose on-shell (OS) renormalisation conditions for the gauge and Higgs sectors, i.e. the wave function renormalisations are chosen to yield unit residues for the propagators at the physical masses M 2 W , M 2 Z , M 2 H ; we remove Z − γ mixing for on-shell kinematics, Σ T AZ (q 2 = 0) = Σ T ZA (q 2 = M 2 Z ) = 0 (where Σ T denotes the transverse part of the self-energy in case of the vector bosons). Alongside on-shell renormalisation in the fermion-sector f (the form of which is not relevant for our work), this reduces the renormalisation condition for the electromagnetic coupling constant to the condition of the amputated three-point function Γ µ ff γ = −ieγ µ for on-shell legs and zero momentum transfer. The renormalisation constants can then be worked out as functions of the self-energies and their derivatives and are tabled in, e.g., [67,72]. From this follows in particular  [75][76][77][78] alongside PackageX [79] for numerical evaluations and analytical cross-checks. Calculating in general gauges ∼ ξ A , ξ W , ξ Z , we have verified gauge independence at the level of divergent parts as well as finite parts (see below). Further checks include comparisons against existing results: In particular Ref. [61] provides a comprehensive overview of the HEFT Lagrangian at one-loop. Given that we find agreement with the results quoted there, we limit ourselves to a summary of the findings relevant for our analysis in Sec. 4.

H → W W * amplitude renormalisation
The (off-shell) H(q)W µ (k 1 )W ν (k 2 ) one-loop three-point function is pictorially represented by the Feynman diagrams of Fig. 2. We limit ourselves to transverse polarisations. The counter term vertex (in a convenient phase convention) follows from the procedure described in Sec. 3.2, see also appendix A, H The remaining singularity ∼ q 2 HW µ W µ then allows us to read off (3.24) The precise form of these additional terms is not relevant for the analysis of Sec. 4, but we quote these here to highlight the agreement with [61]. The singularity related to the bare momentum-independent HW W vertex, on the other hand, fixes

25) after entering the renormalisation constants of appendix A.
H → ZZ * amplitude renormalisation Similar to the discussion above, the (off-shell) H(q)Z(k 1 )Z(k 2 ) one-loop three-point function is summarised in Fig. 3. The counter term vertex is H Z Z f f f Figure 3.  Fig. 3). Again the terms δ[. . . ] enforce non-trivial identities among the one-loop counter terms of the operators of Eq. (3.12), which are again not relevant for our κ 2V − κ V analysis. Entering the wave-function renormalisation constant of appendix A, we obtain a UV-finite on-shell H → γZ amplitude which is manifestly gauge-independent. Our result agrees with the findings of [60] and the unitary gauge calculation of [80].

H → γγ amplitude renormalisation
The H → γγ amplitude has been first computed in general R ξ gauge in [60] after unitary gauge results have been used in [80]. The H → γγ is free of singularities (see also [81,82]) and therefore only imposes one-loop counter term relations among the coefficients of Eq. (3.12). Again these are not relevant for the analysis below, but we have cross-checked our findings against the existing results in the literature.

Decay widths
The amplitudes H → γγ, γZ of the previous section can be converted into decay widths Γ(H → γγ), Γ(H → γZ) straightforwardly; they are only sensitive to ζ 1 . For the cases of H → ZZ * and H → W W * we follow the SMEFT analysis of Refs. [83,84]. To this end, we write the decay width and both the leading order (LO) and one-loop parts can be written as where the polarisation vector of the on-shell boson is denoted with µ . M ρ f f V is the f f V vertex arising from V -current and includes the on-shell spinors, while ∆ νρ V V is the vector boson propagator. We do not include higher order corrections for these parts and only consider them for the HV V part of the amplitude, M µν HV V . In the case of H → ZZ * we also do not include the H → γZ contribution (similar to [84]).
The renormalised HV V vertices described above can be decomposed into three different Lorentz structures using the metric tensor g µν and the Levi-Civita antisymmetric tensor ε µνρσ as The last two terms only contribute at one-loop order for the SM and HEFT for the choices detailed in Sec. 3.2. F NLO 3 is zero for the H → ff Z channel, while the longitudinal-like second term does not contribute as we consider massless fermions. The non-trivial form factors are directly related to the one-loop calculation detailed above. By considering the unpolarised amplitude and performing the integrations of Eq. (3.29) we verify that our calculations reproduce the expressions supplied by Refs. [83][84][85] for the SM case at LO (see also [86,87]). For HEFT with NLO corrections, the integrations are performed numerically to obtain the decay widths Γ(H → W W * ) and Γ(H → ZZ * ) and subsequently the relevant κ V = 1 + ζ 1 for a scale choice µ R = M H . For an efficient evaluation of the amplitudes, we furthermore choose Feynman Gauge ξ W = ξ Z = ξ A = 1 and rescale our results to the values reported by the Higgs Cross Section Working Group [88] for the SM point ζ 1 = ζ 2 = 0 for ease of comparison.

Comments on κ 2V − κ V interpretations
The cancellation of gauge fixing parameters ξ W,Z,A in the amplitude calculations detailed above is highly non-trivial. Indeed, individual Feynman diagrams of Figs. 2 and 3 contain singularities ∼ ξ 2 W,Z , ξ W,Z,A , which cancel when we sum over the contributing diagrams. For instance, the vanishing ξ 2 W,Z ∆ UV dependence of H(q) → W (k 1 )W (k 2 ) results from a cancellation between the gauge boson insertions of diagrams Fig. 2 (g) and (h) like in the SM. A zero ∼ ξ W,Z ∆ UV coefficient results from cancellations between diagrams Fig. 2 (b), (c), (d), (i) and (j), and ∼ ξ A ∆ UV is due to destructive interference between (c), (g), (j). Similar cancellations take place for H → ZZ (there is dependence on ξ A ).
As can be seen from Eqs. (3.23) and (3.24) (see also appendix A), terms of Eq. (3.12) are sourced even for ζ 1 = ζ 2 = 0. Although we must include all relevant HEFT coefficient renormalisation factors in the calculation, we can obtain a consistent κ 2V − κ V correlation at our scale choice µ R = M H by projecting all amplitudes to a i = 0 after renormalisation. Concretely this means that we will limit our attention to ζ 1,2 in the following and choose all other a i (M H ) = 0 as MS parameters. In particular, this means that operator matrix elements HV µν V µν (V = A, W, Z) are absent, which would otherwise impact the Higgs boson production and decay phenomenology. These should be included in a comprehensive fit of the HEFT Lagrangian that we do not attempt here.

κ V − κ 2V correlations from Higgs data
Equipped with the findings of the previous sections, we are now ready to turn to results. For the current status of the Higgs programme, we construct a χ 2 test from the κ results provided in Ref. [89]. The χ 2 statistic defined from the κ data at an integrated luminosity of 139 fb −1 is given by where κ i,exp are the central values and V ij is the covariance matrix, the elements of which are given as ρ ij σ i σ j with σ i as uncertainties and ρ ij as correlations. The details of the experimental data are provided in Tab. 2. Using the four κ parameters detailed above in Eq. (4.1), we perform a χ 2 fit for ζ 1 and ζ 2 . 6 The allowed regions are shown in Fig. 4a.
To gauge the improvements that can be expected at the high luminosity (HL) LHC phase, we extrapolate these contours to an integrated luminosity of 3000 fb −1 . To achieve this, we scale the ATLAS uncertainties for the current luminosity (Lumi ATLAS ) of 139 fb −1 by a naive scaling factor The projected HL-LHC uncertainties are tabulated in column 3 of Tab. 2 for an unchanged correlation matrix (column 4). We stress that we use this extrapolation purely to obtain a qualitative outlook for the HL-LHC, but note that this rescaling provides a good approximation of the comprehensive investigation of [90]. The resulting constraints at 68% and 95% confidence, assuming consistency with the SM, are shown in Fig. 4b. As is visible from Fig. 4, we have rescaled the ζ 2 axes to account for the loop suppression. Given that the correlation κ V (κ 2V ) is related to a weak radiative correction, the bounds on κ 2V = 1 + ζ 2 from single Higgs data are loose, even when we consider the direct sensitivity of the early stage LHC programme shown in Fig. 1. Extrapolations to 3/ab show that single Higgs physics will greatly collapse the current limits along the ζ 1 direction, however the loop suppression implies a sensitivity |ζ 2 | ∼ 4.  Fig. 4, based on [89]. Columns 2 and 4 show the ATLAS Run-2 results along with the correlation matrix. In column 3, the HL-LHC projected uncertainties are listed.  Figure 4. The ζ 1 vs. ζ 2 parameter space obtained from the χ 2 analysis of the κ parameters from ATLAS Run-2 data [89]. The solid (dashed) contour denotes the ∆χ 2 = 2.28 (∆χ 2 = 5.99) which corresponds to confidence level of 68% (95%). Fig. 4b shows the similar regions in blue obtained with the HL-LHC projected data.
A coupling modifier of this size is already in tension with the current measurement results detailed in Sec. 2. This shows that the κ 2V and κ V are largely independent parameters at the LHC as is naively suggested by Eq. (3.5). This highlights the need to increase the sensitivity coverage to ζ 2 through direct searches further, which we will turn to in Sec. 5.

Enhancing direct collider sensitivity with Graph Neural Networks
Higgs pair production via WBF has been discussed in Refs. [21][22][23][24][25] (see also [91] for a recent dimension 8 SMEFT study). The process shares the properties of single Higgs WBF production [92][93][94][95][96] especially when we consider LHC-relevant QCD corrections [97][98][99]. To select the events for the WBF topology, Fig. 5, we choose two forwarded jets in opposite detector hemispheres (opposite signs of pseudorapidity η j1 · η j2 < 1) with the invariant mass m jj larger than 500 GeV. For the minimum transverse momentum p T of WBF jets, we Figure 5. Representative Feynman diagrams contributing to WBF p(f )p(f ) → HHjj. We consider H → bb final states, but do not include the decays to the diagrams.  Figure 6. A representative diagram of an event as a fully connected graph is shown on the left. The network performance is shown on the right with a receiver operating characteristic (ROC) curve along with the corresponding AUC=0.995. choose 50 GeV. We select events with four b-jets in the central region |η b−jets | < 2.5 with minimum p T ≥ 20 GeV. 8 Events are generated in proton-proton collisions at √ s = 13 TeV using MadGraph5 aMC@NLO [100] at leading order precision using the standard setting scale choices. 9 The main SM background for the 4b + 2j final state comes from multijet QCD processes [21][22][23][24][25]. After the given VBF selection cuts, the QCD multijet background cross-section is 4.41 pb, while the signal cross-section is 0.086 f b, which highlights the necessity to efficiently reduce the background.
Going beyond 'traditional' techniques, utilising Graph Neural Networks (GNNs) is a motivated avenue. These strategies are tailored to accessing the graph structure of differential cross section data under realistic conditions that we typically parametrise through Feynman diagram calculations in theoretical reference calculations. Hence, recent applications to particle physics data in a range of areas stretching from QCD [101], over anomaly detection [102] (including IR safe formulations [103,104]), to SMEFT parameter fits [105] have highlighted GNNs as versatile and highly sensitive tools for BSM discrimination. Given the special phenomenological properties of WBF processes, we can therefore expect similar performance improvements in this channel.
To construct a GNN for the present signal vs. background discrimination task, we use a fully-connected bi-directional graph to represent the 4b + 2j event (Fig. 6a). Each jet is 8 We use a flat tagging efficiency of 70% and ignore mistagging. For the EdgeConv GNN strategy we employ below, mistagging will predominantly change the background normalisation, which can be compensated through a tuned working point of the classifier. We comment on uncertainties from the background normalisation further below. 9 It is known that the NLO QCD corrections are relatively mild (e.g. [93,98]) and can be largely captured through adapted weak-boson fusion factorisation scale choices inspired by the "double-DIS" structure of WBF [94,95]. We will see below that the used approach reproduces the current experimental outcome reasonably well to warrant a high-luminosity extrapolation.
represented by a node, which is associated with a node feature vector [p T , η, φ, E, m, P ID] (representing transverse momentum, pseudorapidity, azimuthal angle, energy, mass and particle identification number, respectively). To implement the graph structure, we use the Deep Graph Library [106] and PyTorch [107] and choose an Edge Convolution (EdgeConv) network to classify signal and background. Edge convolution is known to be particularly suited for extracting high-level variable information as edge features from given low-level node features [108]. The message passing function for the Edge convolution is defined as x Here, x (l) i represents the input node features of node i in the l-th message passing layer, with l = 0 denoting the input node features of the graph. The neighbourhood set N (i) consists of all nodes in the graph connected to node i. The linear layers Θ and Φ take the vectors x j − x i and x i , respectively, and map them to the same vector space. We apply a single message passing step producing 40-dimensional node features. 10 Since we aim to classify graphs, we apply a mean graph readout operation to these node features. A classifier multilayer perceptron with single layers having 40 nodes, with ReLU activation, takes these features and outputs a vector of two dimensions (after being properly normalised with a SoftMax activation function). We minimise the cross-entropy loss function using the Adam optimiser [109] with a learning rate of 10 −3 . Given the distinct phenomenological properties, we also find that already a shallow network is well-suited to obtain good discrimination.
The network is trained for 100 epochs with batch sizes of 124 for a κ 2V -enriched parameter point sample (κ 2V , κ V = 2, 1). We split the total data into 80%, 10% and 10% for training, validation and testing, respectively. We verify that the validation and the training accuracies are comparable to avoid overtraining. The receiver operating characteristic curve of the network is shown in Fig. 6b with an area under curve (AUC) of 0.995, implying a very high classifier performance. 11 With the classifier optimised as described above, we can obtain the limit on κ 2V by fixing κ V = 1 (we do not include modified trilinear Higgs interactions as in the previous sections) to compare the performance of the GNN with ATLAS analysis of Ref. [37]. For the optimal working point on the ROC curve, the QCD multijet background cross-section drops to is 0.88 fb, while the signal cross-section remains 0.037 fb. The exclusion limits are calculated using the significance σ = N s / √ N s + N b , where N s is the number of signal events and N b is the number of background events at a particular luminosity. This is shown in Fig. 7a, where we show the ATLAS constraint as a horizontal (cyan) line. Firstly, this agreement acts as a good validation of our procedure in comparison to the realistic experimental measurement, which gives us confidence for the HL extrapolation that we will perform next. The relative improvement suggested in Fig. 7a might be optimistic. Given  Figure 7. Exclusion limits on κ 2V and κ V parameter space obtained from the GNN output trained for the κ 2V enriched sample (κ 2V , κ V = 2, 1) are shown in Fig. 7a. The cyan line represents the limits obtained from the ATLAS analysis for the fixed value of κ V = 1 [37]. The dot-dashed contour shows the effect of increasing the background by 25%. In Fig. 7b, HL-LHC projected constraints are overlaid on the 95% C.L. region shown in Fig. 4b for κ 2V and κ V parameter space.
the necessity to overcome large backgrounds, the choice of working points and the accuracy of the background simulation are factors that can impact the discrimination, in particular when we deal with the large statistics of the HL phase. The competitiveness of our idealised GNN analysis, however, demonstrates that such techniques deserve consideration as part of realistic experimental analyses. We will also comment on a range of potential architecture improvements that we have not included in our analysis below. Using the κ 2V = 2 as a reference point also for a more general scan, we choose a fixed working point of the classifier which maximises the significance to obtain exclusion contours in the κ 2V , κ V parameter space also shown in Fig. 7a (this includes the changes to the H → bb branching ratio as a function of κ V ). To compare the level of sensitivity, we overlay these, projected to 3/ab HL-LHC target with the expected sensitivity from single Higgs measurements detailed in Sec. 4 in Fig. 7b. Note that in this exploratory comparison, we have not included subdominant background from weak multi-boson processes or top production, which will additionally degrade the sensitivity (see e.g. [22]). These processes are, however, also characterised by different kinematic and QCD radiation properties, and we can expect significant discrimination when considering these additional contributions in the GNN classification. To get an estimate of how additional backgrounds modify the direct search sensitivity, we modify the considered dominant QCD through a naive shifting by 25% in normalisation. Including these to Fig. 7a, we see that there is minor change in the sensitivity, and good sensitivity to κ 2V can be achieved.
The GNN could be further improved through a more dedicated inclusion of κ 2V − κ V correlations to the classification, e.g., through a multi-class network architecture discussed in Ref. [105], which could also be extended to subdominant background contributions to specifically combat those when the QCD contribution has been removed sufficiently. A more robust graph embedding could help in such an approach as well. While the QCD background phenomenology is very different from the WBF signal (highlighted by the fact that a shallow network is sufficient for very good discrimination), different electroweak correlation structures such as top and weak boson decays impart a more tree-like structure that could be further exploited when these become relevant. Insensitivity to differential distributions of the node features, as well as overall normalisations, could be achievable using adversarial networks [110]. We leave more detailed investigations for future work.

Summary and Conclusions
In this work, we have presented a comprehensive and theoretically consistent discussion of κ V − κ 2V correlations by employing Higgs Effective Field Theory, which puts related analyses at the LHC [36-38] on a theoretically firm footing. We show that single Higgs constraints can be formulated for κ 2V analyses, however, given that these effects arise as a weak radiative correction, the single-Higgs κ V constraints are relatively loose and not competitive when compared to the LHC's sensitivity to κ 2V already at this stage of the programme. Indeed, the LHC's sensitivity pattern when mapped onto the HEFT Lagrangian Eq. (3.5) and its weak corrections allows us to treat different Higgs interactions in the gauge sector as largely independent parameters. As κ 2V analyses are mainly driven by the direct investigation of weak HHjj production, an enhanced direct sensitivity of future LHC runs provides the best motivated avenue to obtain a more fine-grained picture of the Higgs boson's gauge interactions along these lines. To this end, we employ Graph Neural Network techniques to demonstrate that the direct sensitivity to κ 2V can be enhanced as GNNs formidably exploit the particular structure of the decay, colour and kinematical correlations of the HHjj signal compared to the dominant QCD backgrounds. Improvements to κ 2V ∼ 1 ± 40% could then be within the reach of the HL-LHC, in particular as improvements through multi-class GNNs and additional architecture improvements are likely to add further sensitivity.

A UV divergent parts of renormalisation constants
In this section, we provide the explicit form of the UV divergent parts renormalisation constants used to renormalise the Lagrangian of Sec. 3.1 required in this work. The field strength renormalisations are obtained by evaluating all relevant (off-shell, one-particle irreducible) two-point function X → Y contributions derived from Eq. (3.5), which we represent as to be contrasted with 2-point insertions from Eq. (3.12) represented by The counter terms are then pictorially determined from the renormalisation conditions (RCs) through The counter term diagrams will typically contain terms from Eq. (3.12) which we highlight via In the conventions detailed above, the contribution to the off-shell Higgs boson two-point function H(q) → H(q) that needs to be included is  Note that the Higgs renormalisation is manifestly gauge-independent as the Higgs is a singlet field in HEFT [61] (this also applies to the tadpole). The transverse part of the W (q) → W (q) boson polarisation function receives no contributions from Eq. ( with