Pseudo-observables in electroweak Higgs production

We discuss how the leading electroweak Higgs production processes at the LHC, namely vector-boson fusion and Higgs+W/Z\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$+W/Z$$\end{document} associated production, can be characterized in generic extensions of the Standard Model by a proper set of pseudo-observables (PO). We analyze the symmetry properties of these PO and their relation with the PO set appearing in Higgs decays. We discuss in detail the kinematical studies necessary to extract the production PO from data, and present a first estimate of the LHC sensitivity on these observables in the high-luminosity phase. The impact of QCD corrections and the kinematical studies necessary to test the validity of the momentum expansion at the basis of the PO decomposition are also discussed.


Introduction
Characterizing the properties of the Higgs boson, both in production and in decay processes, with high precision and minimum theoretical bias, is one of the main goal of future experimental efforts in high-energy physics and a promising avenue to shed light on physics beyond the Standard Model (SM). In this context, a useful tool is provided by the so-called Higgs pseudo-observables (PO) [1][2][3][4][5]. The latter constitute a finite set of parameters that are experimentally accessible, are well defined from the point of view of quantum field theory (QFT), and characterize possible deviations from the SM in processes involving the Higgs boson in great generality. More precisely, the Higgs PO are defined from a general decomposition of on-shell amplitudes involving the Higgs boson -based on analyticity, unitarity, and crossing symmetry -and a momentum expansion following from the dynamical assumption of no new light particles (hence no unknown physical poles in the amplitudes) in the kinematical regime where the decomposition is assumed to be valid.
The idea of PO has been formalized the first time in the context of electroweak observables around the Z pole [6,7], while the generalization relevant to analyze Higgs decays has been presented in Ref. [1]. In this paper we further generalize the PO approach to describe electroweak Higgs-production processes, namely vector-boson fusion (VBF) and associated production with a massive SM gauge boson (VH).
The interest of such production processes is twofold. On the one hand, they are closely connected to the h → 4 , 2 2ν decay processes by crossing symmetry, and by the exchange of lepton currents into quark currents. As a result, some of the Higgs PO necessary to describe the h → 4 , 2 2ν decay kinematics appear also in the description of the VBF and VH cross sections (independently of the Higgs decay mode). This fact opens the possibility of combined analyses of production cross sections and differential decay distributions, with a significant reduction on the experimental error on the extraction of the PO. On the other hand, studying the production cross sections allows us to explore different kinematical regimes compared to the decays. By construction, the momentum transfer appearing in the Higgs decay amplitudes is limited by the Higgs mass, while such a limitation is not present in the production amplitudes. This fact allows us to test the momentum expansion that is intrinsic in the PO decomposition, as well as in any effective field theory approach to physics beyond the SM.
Despite the similarities at the fundamental level, the phenomenological description of VBF and VH in terms of PO is significantly more challenging compared to that of Higgs decays. On the one hand, QCD corrections play a nonnegligible role in the production processes. Although technically challenging, this fact does not represent a conceptual problem for the PO approach: the leading QCD corrections factorize in VBF and VH, similarly to the factorization of QED corrections in h → 4 [8]. As we will show, this implies that NLO QCD corrections can be incorporated in general terms with suitable modifications of the existing Monte Carlo tools. On the other hand, the relation between the kinematical variables at the basis of the PO decomposition (i.e. the momentum transfer of the partonic currents, q 2 ) and the kinematical variables accessible in pp collisions is not straightforward, especially in the VBF case. As we will show, this problem finds a natural solution in the VBF case due to strong correlation between q 2 and the p T of the VBF-tagged jets.
The paper is organized as follows: in Sect. 2 we present the decomposition in terms of PO of the electroweak amplitudes relevant to VBF and VH, analyzing the relation with the decay PO already introduced in Ref. [1]. In Sect. 3 we present a phenomenological analysis of the VBF process, discussing in detail the implementation of QCD corrections, and the key role of the jet p T for the identification of the PO. An estimate of the statistical error expected on the PO extracted from VBF in the high-luminosity phase at the LHC is also presented. A similar discussion for the VH processes is presented in Sect. 4. A detailed discussion as regards the validity of the momentum expansion, and how to test it from data, is presented in Sect. 5. The results are summarized in Sect. 6.

Amplitude decomposition
Neglecting light fermion masses, the electroweak production processes VH and VBF or, more precisely, the electroweak partonic amplitudes f 1 f 2 → h + f 3 f 4 , can be completely described by the three-point correlation function of the Higgs boson and two (color-less) fermion currents where all the states involved are on-shell. The same correlation function also controls four-fermion Higgs decays.
In the h → 4 , 2 2ν case both currents are leptonic and all fermions are in the final state [1]. In case of VH associate production one of the currents describes the initialstate quarks, while the other describes the decay products of the (nearly on-shell) vector boson. Finally, in VBF production the currents are not in the s-channel as in the previous cases, but in the t-channel. Strictly speaking, in VH and VBF the quark states are not on-shell; however, their offshellness can be neglected compared to the electroweak scale characterizing the hard process (both within and beyond the SM). Following Ref. [1], we expand the correlation function in Eq. (1) around the known physical poles due to the propagation of intermediate SM electroweak gauge bosons. The PO are then defined by the residues on the poles and by the nonresonant terms in this expansion. By construction, terms corresponding to a double pole structure are independent from the nature of the fermion current involved. As a result, the corresponding PO are universal and can be extracted from any of the processes mentioned above, both in production and in decays.

Vector-boson fusion Higgs production
Higgs production via vector-boson fusion (VBF) receives contribution both from neutral-and charged-current channels. Also, depending on the specific partonic process, there might be two different ways to construct the two currents, and these two terms interfere with each other. For example, in uu → uuh two neutral-current processes interfere, while in ud → udh there is an interference between neutral and charged currents. In this case it is clear that one should sum the two amplitudes with the proper symmetrization, as done in the case of h → 4e [1].
We now proceed describing how each of these amplitudes can be parametrized in terms of PO. Let us start with the neutral-current one. The amplitude for the on-shell process where q 1 = p 1 − p 3 , q 2 = p 2 − p 4 , and T μν n.c. (q 1 , q 2 ) is the same tensor structure appearing in h → 4 f decays [1]. In particular, Lorentz invariance allows for only three possible tensor structures, to each of which we can assign a generic form factor: The form factor F L describes the interaction with the longitudinal part of the current, as in the SM; the F T term describes the interaction with the transverse part, while F C P describes the CP-violating part of the interaction (if the Higgs is assumed to be a CP-even state). The charged-current contribution to the amplitude for the on-shell process where, again, T μν c.c. (q 1 , q 2 ) is the same tensor structure appearing in the charged-current h → 4 f decays: The amplitudes for the processes with anti-quarks in the initial state can easily be obtained from the above ones. The next step in the decomposition of the amplitude requires one to perform a momentum expansion of the form factors around the physical poles due to the propagation of SM electroweak gauge bosons (γ , Z , and W ± ), and to define the PO (i.e. the set {κ i , i }) from the residues of such poles. We stop this expansion neglecting terms which can be generated only by local operators with dimension higher than six. A discussion as regards limitations and consistency checks of this procedure is presented in Sect. 5. The explicit form of the expansion of all the form factors in terms of PO can be found in Ref.
[1] 1 and will not be repeated here. We report here explicitly only expressions for the longitudinal form factors, which are the only ones containing PO not present also in the leptonic decay amplitudes: 1 With respect to [1] we modified the labels of the form factors: F 1 → F L , F 3 → F T , and F 4 → F C P , and analogously for the G i . Here is the sine (cosine) of the Weinberg angle. 2 The functions SM L ,n.c.(c.c.) (q 2 1 , q 2 2 ) denote non-local contributions generated at the one-loop level (and encoding multi-particle cuts) that cannot be re-absorbed into the definition of κ i and i . At the level of precision we are working, taking into account also the high-luminosity phase of the LHC, these contributions can safely be fixed to their SM values.
As anticipated, the crossing symmetry between h → 4 f and 2 f → h 2 f amplitudes ensures that the PO are the same in production and decay (if the same fermions species are involved). The amplitudes are explored in different kinematical regimes in the two type of processes (in particular the momentum transfer, q 2 1,2 , are space-like in VBF and timelike in h → 4 f ). However, this does not affect the definition of the PO. This implies that the fermion-independent PO associated to a double pole structure, such as κ Z Z and κ W W in Eq. (6), are expected to be measured with higher accuracy in h → 4 and h → 2 2ν rather than in VBF. On the contrary, VBF is particularly useful to constrain the fermiondependent contact terms Zq i and W u i d j , which appear only in the longitudinal form factors. For this reason, in the following phenomenological analysis we focus our attention mainly on the LHC reach on these parameters. Still, we stress that the PO framework is well suited to perform a global fit including production and decay observables at the same time.

Associated vector-boson plus Higgs production
By VH we denote the production of the Higgs boson with a nearly on-shell massive vector boson (W or Z ), starting from and initial qq state. For simplicity, in the following we will assume that the vector boson is on-shell and that 2 More precisely, (g ik the interference with the VBF amplitude can be neglected. However, we stress that the PO formalism clearly allows one to describe both these effects (off-shell V and interference with VBF in the case of V →qq decay) simply applying the general decomposition of neutral-and charged-current amplitudes as outlined above. Similarly to VBF, Lorentz invariance allows us to decompose the amplitudes for the on-shell processes in three possible tensor structures: a longitudinal one, a transverse one, and a CP-odd one, where q = p 1 + p 2 = k + p. In the limit where we neglect the off-shellness of the final-state V , the form factors can only depend on q 2 . Already from this decomposition of the amplitude it should be clear that differential measurements of the VH cross sections as a function of q 2 [9], as well as in terms of angular variables that allow one to disentangle different tensor structures, are an important input to constrain the PO. Performing the momentum expansion of the form factors around the physical poles, and defining the PO as in Higgs decays and VBF, we find where we have omitted the indication of the (tiny) non-local terms, fixed to their corresponding SM values. According to the arguments already discussed at the end of Sect. 2.1, in the following phenomenological analysis we focus our attention on the longitudinal form factors F L and G L and, in particular, on the extraction of the quark contact terms Zq i and W u i d j .
2.3 Parameter counting, symmetry limits, and dynamical assumptions on the PO We now want to analyze the number of free parameters and the symmetry limits for the newly introduced PO appearing in VBF and VH production, compared to the decay PO introduced in Ref. [1]. The additional set of PO (the "production PO") is represented by the contact terms for the light quarks.
In a four-flavor scheme, in absence of any symmetry assumption, the number of independent parameters for the neutralcurrent contact terms is 16  denotes the two quark doublets in the basis where down quarks are diagonal. This symmetry is an exact symmetry of the SM in the limit where we neglect light quark masses. Enforcing it at the PO level is equivalent to neglecting terms that do not interfere with SM amplitudes in the limit of vanishing light quark masses. Under this (rather conservative) assumption, the number of independent neutral-current contact terms reduces to eight real parameters, 3 (10) and only two complex parameters in the charged-current case: A further interesting reduction of the number of parameters occurs under the assumption of an U (2) 3 symmetry acting on the first two generations, namely the maximal flavor symmetry compatible with the SM gauge group [12][13][14].
The independent parameters in this case reduces to six: where W u L is complex, or five if we further neglect CPviolating contributions (in such case W u L is real). We employ this set of assumptions (U (2) 3 flavor symmetry and CP conservation) in the phenomenological analysis of VBF and VH processes discussed in the rest of the paper. Finally, we can enforce custodial symmetry that, as shown in [1], implies reducing the number of independent PO to four in the U (2) 3 case (independently of any assumption as regards CP). As far as dynamical hypotheses are concerned, numerical constraints on the Higgs PO can be derived under the hypothesis that the Higgs particle is the massive excitation of a pure SU (2) L doublet, i.e. within the so-called linear EFT (or SMEFT). In this framework the Higgs PO receive contributions from effective operators written in terms of the doublet field H , which contribute also to non-Higgs observables. As a result, it is possible to derive relations between the Higgs PO and electroweak precision observables, as well as relations among Higgs PO which reduce the number of independent parameters. The matching to the SMEFT at the dimension-6 level in various bases, and the explicit relations among Higgs PO which follow, can be found in Refs. [1,2]. Limiting the attention to the (presumably dominant) treelevel contributions, generated by dimension-6 operators, the following relations can be derived [2]: where δg Z f and δg W f are the effective Z -and W -couplings to SM fermions, δg 1,z and δκ γ are the anomalous triple gauge couplings (aTGC), and T 3 f and Y f are the isospin and hypercharge quantum numbers of the fermion f . Moreover, the custodial-symmetry relation (13) is automatically enforced at the dimension-6 level.
Recent analyses of Z -and W -pole observables within the SMEFT, with a generic flavor structure, can be found in Refs. [15,16]. A combined fit to LEP-II WW and LHC Higgs signal strengths data, which removes all the flat directions in the determination of aTGC within the SMEFT has been presented in Ref. [17]. Combining some of these recent fits (in particular Z -and W -pole couplings from Ref. [15] and aTGC from Ref. [17]) we find the following numerical constraints on the quark contact terms (within the SMEFT): where, for simplicity, we have further imposed the U (2) 3 flavor symmetry hypothesis. The corresponding correlation matrix turns out to be close to the identity matrix. The precise values of these results is not relevant to the present analysis, but it can be used as a guideline for the sensitivity needed on the PO measured from VBF and VH in order to tests SMEFT predictions. As we show later on, the LHC at high luminosity will reach such a sensitivity.
A further restrictive dynamical hypothesis is obtained within the framework of the so-called universal theories, i.e. by assuming that all new physics interactions can be written in terms of the SM bosonic fields only. All analyses of VBF and VH production, as well as of h → 4 f decays, performed assuming new physics only via modified hV V vertices belong to this category, e.g. in Refs. [18][19][20][21]. A specific example of this scenario are the parametric expressions of the Higgs PO in terms of the socalled "Higgs characterization framework" introduced in Refs. [18,19] 4 : In this case the variability of the neutral-current contact terms is further reduced by a dynamical assumption that links them to the two terms κ H ∂ Z and κ H ∂γ . We stress that such an assumption cannot be justified only in terms of symmetry principles.
Using FeynRules [10] we implemented a general UFO model [11] containing all the Higgs PO (including also decays [1]). The model itself will promptly be made available online [22] and allows for comprehensive phenomenological Monte Carlo studies at the LHC. A detailed implementation of the Higgs PO framework in a Monte Carlo tool including NLO QCD corrections will be presented in a subsequent publication.

VBF kinematics
Vector-boson fusion Higgs production is the largest of all electroweak Higgs-production mechanisms in the SM at the LHC. It is highly relevant in the context of experimental Higgs searches due to its striking signature, i.e. two highly energetic forward jets in opposite detector hemispheres, which allows an effective separation from the backgrounds. In this chapter we study the phenomenology of VBF production in the PO framework. We mainly concentrate our discussion on measuring the quark contact term PO, Zq i and W u i d j , namely the residues of the single pole terms in the expansion of the longitudinal form factors in Eq. (6).
At the parton level (i.e. in the qq → hqq hard scattering) the ideal observable relevant to extract the momentum dependence of the factor factors would be the double differential cross section d 2 σ/dq 2 1 dq 2 2 , where q 1 = p 1 − p 3 and q 2 = p 2 − p 4 are the momenta of the two fermion currents entering the process (here p 1 , p 2 ( p 3 , p 4 ) are the momenta of the initial (final) state quarks). These q 2 i are also the key variables to test and control the momentum expansion at the basis of the PO decomposition.
As a first step of the VBF analysis we have to choose a proper pairing of the incoming and outgoing quarks, given we are experimentally blind to their flavor. For partonic processes receiving two interfering contributions when the final-state quarks are exchanged, such as uu → huu or ud → hud, the definition of q 1,2 is even less transparent since a univocal pairing of the momenta cannot be assigned, in general, even if one knew the flavor of all partons. This problem can be overcome at a practical level by making use of the VBF kinematics, in particular by the fact that the two jets are always very forward. This implies that one can always pair the momentum of the jet going in one direction with the initial parton going in the same direction. In the same way we can argue that the interference between different amplitudes (e.g. neutral current and charged current) is negligible in VBF. In order to check this, we perform a leading order (LO) parton level simulation of VBF Higgs production ( pp → h j j at O(α 3 )) employing MadGraph5_aMC@NLO [23] (version 2.2.3) at 13 TeV c.m. energy together with the Higgs PO UFO model. In this simulation we impose the basic set of VBF cuts, p T,j 1,2 > 30 GeV, |η j 1,2 | < 4.5, and m j 1 j 2 > 500 GeV.
In Fig. 1, we show the distribution in the opening angle of the incoming and outgoing quark momenta for the two different pairings. The left plot shows the SM, while the right plot shows a specific NP benchmark point. Depicted in blue is the pairing based on the leading color connection using the color flow variable in the event file, while in red we show the opposite pairing. The plot shows that the momenta of the color connected quarks tend to form a small opening angle and the overlap between the two curves, i.e. where the interference effects might be sizable, is negligible. This implies that in the experimental analysis the pairing should be done based on this variable. Importantly, the same conclusions can be drawn in the presence of new physics contributions to the contact terms.
There is a potential caveat to the above argument: the color flow approximation ignores the interference terms that are of higher order in 1/N C . Let us consider a process with two interfering amplitudes with the final-state quarks exchanged, for example in uu → uuh. The differential cross section receives three contributions propor- For the validity of the momentum expansion it is important that the momentum transfers (t i j ) remain smaller than the hypothesized scale of new physics. On the other hand, imposing the VBF cuts, the interference terms turn out Fig. 2 Leading-order parton level simulation of the Higgs VBF production at 13 TeV pp c.m. energy. Shown here is the density histogram in two variables; the outgoing quark p T and the momentum transfer −q 2 with the initial "color-connected" quark. The left plot is for the SM, while the plot on the right is for a specific NP benchmark to depend on one small and one large momentum transfer. However, thanks to the pole structure of the form factors, they give a very small contribution.
Even though in some experimental analyses, after reconstructing the momenta of the two VBF-tagged jets and the Higgs boson, one could in principle compute the relevant momentum transfers q 1 and q 2 , adopting the pairing based on the opening angle, in an hadron collider environment like the LHC this is unfeasible. Furthermore, for other Higgs decays modes, such as h → 2 2ν, it is not possible to reconstruct the Higgs boson momentum. Therefore, we want to advocate the use of the p T of the VBF jets as a proxy for the momentum transfers q 2 1,2 . The quality of this approximation can be understood by explicitly computing the momentum transfers q 2 1,2 in the VBF limit | p T | E jet and for a Higgs produced close to threshold. Let us consider the partonic momenta in the c.o.m. frame for the process: p 1 = (E, 0, E), . Conservation of energy for the whole process dictates 2E = E 1 + E 2 + E h , where E h is the Higgs energy, usually of order m h if the Higgs is not strongly boosted. In this case E − E i = E i E since the process is symmetric in 1 ↔ 2. For each leg, energy and momentum conservation (along the z axis) give Putting together these two relations, one finds where in the last step we assumed E i E , i.e. the Higgs being produced near threshold.
In order to confirm the above conclusion, in Fig. 2 we show a density histogram in two variables: the (observable) p T of the outgoing jet and the (unobservable) momentum transfer −q 2 obtained from the correct color flow pairing (the left and the right plots are for the SM and for a specific NP benchmark, respectively). These plots indicate a very strong correlation of the jet p T with the momentum transfer −q 2 associated with the correct color pairing. We stress that this conclusion holds both within and beyond the SM. Therefore, we encourage the experimental collaborations to report the unfolded measurement of the double differential distributions in the two VBF-tagged jet p T :F( p T j 1 , p T j 2 ). This measurable distribution is indeed closely related to the form factor entering the amplitude decomposition, F L (q 2 1 , q 2 2 ), and encodes (in a model-independent way) the dynamical information as regards the high-energy behavior of the process. Moreover, as we will discuss in Sect. 3.3, the extraction of the PO in VBF must be done preserving the validity of the momentum expansion: the latter can be checked and enforced setting appropriate upper cuts on the p T distribution. As an example of the strong sensitivity of the (normalized)F( p T j 1 , p T j 2 ) distribution to NP effects, in Fig. 3, we show the corresponding prediction in the SM (left plot) and for a specific NP benchmark (right plot).

NLO QCD corrections in VBF
Inclusive VBF Higgs production in the SM is very stable with respect to higher-order QCD corrections [24][25][26][27]. Employing a fixed renormalization and factorization scale μ R,F = m W inclusive NLO QCD corrections are at the level of 5-10 % with remaining scale uncertainties of a few percent. At the NNLO QCD level these uncertainties on the inclusive cross section are further reduced below 1 % [28,29]. However, in more exclusive observables, like the p T spectra of the VBF jets, or when more exclusive experimental selection cuts are applied, sensitivity to QCD radiation is more severe [25], yielding non-negligible NLO correction factors while NLO scale uncertainties remain small (mostly well below 10 %). It is then important to take into account such effects when analyzing VBF beyond the SM, as illustrated in detail in Ref. [19]. Recently the dominant NNLO QCD corrections have been calculated fully differentially [30] pointing towards a non-trivial phase-space dependence with 5-10 % corrections with respect to NLO. Besides higher-order corrections of QCD origin, also EW corrections are relevant for VBF Higgs production [31,32]. At an inclusive level they amount to about −5 % [31], while at the differential level due to the presence of large EW Sudakov logarithms they reach for example −15 % for p T,j 1 = 400 GeV and −10 % for p T,j 2 = 150 GeV [32].
In the following we will illustrate that the perturbative convergence for exclusive VBF observables can be improved when using a dynamical scale μ 0 = H T /2 (with H T being the scalar sum of the p T of all final-state particles) with respect to a fixed scale μ 0 = m W . In particular, here we will focus on the p T spectra of the VBF jets -as inputs for a fit of the Higgs PO. To this end we employ the fully automated Sherpa+OpenLoops framework [33][34][35][36][37][38] for the simulation of EW production of pp → h j j at LO and NLO QCD in the SM. Before applying the VBF selection cuts defined in Eq. (17) we cluster all final-state partons into anti-k T jets with R = 0.4 and additionally require a rapidity separation of the two hardest jets of η j 1 j 2 > 3. This additional requirement, could slightly reduce the capability of differentiating different tensor structures [19], however, such a cut is, on the one hand, experimentally required in order to suppress QCD backgrounds. 5 On the other hand, without such a cut NLO predictions for the p T spectra of the jets become highly unstable when the VBF jet selection is just based on the hard- 5 In fact, in most VBF analyses an even tighter selection of η j1 j2 > 4.5 is imposed. ness of the jets, i.e. a bremsstrahlung jet is easily amongst the two hardest jets and spoils the correlation between the p T of the jets and the momentum transfer, as discussed in Sect. 3.1.
Thanks to the dynamical scale choice NLO corrections to the one-dimensional distributions are almost flat and amount to about −15 %, while the dependence in the twodimensional distribution remains moderate with largest corrections for p T,j 1 ≈ p T,j 2 .
In the following section we will detail a fit of Higgs PO based on LO predictions of VBF using the scale choice and setup developed in this chapter. Here we already note that this fit is hardly affected by the overall normalization of the predictions. Thus, with respect to possible small deviations from the SM due to effective form factor contributions we expect a very limited sensitivity to QCD effects assuming a similar stabilization of higher-order corrections as observed for the SM employing the scale choice μ 0 = H T /2.
In order to verify this assumption and to improve on the Higgs PO fit, we are currently extending the simulations within the Higgs PO framework to the NLO QCD level. To this end, the framework has been implemented in the Open-Loops one-loop amplitude generator in a process independent way. Here, the O(α S ) rational terms of R 2 -type required in the numerical calculation of the one-loop amplitudes in OpenLoops have been obtained generalising the corresponding SM expressions [40]. The implementation of the dipole subtraction and partonshower matching in the Sherpa Monte Carlo framework is based on the model independent UFO interface of Sherpa [41] and is currently being validated.

Prospects for the Higgs PO in VBF at the HL-LHC
The extraction of the PO from the double differential dis-tributionF( p T j 1 , p T j 2 ) has to be done with care. Here we make an attempt to perform such an analysis. In the following we estimate the sensitivity of the HL-LHC, operated at 13 TeV with 3000 fb −1 of data, on measuring the PO assuming maximal flavor symmetry in a seven-dimensional fit to κ Z Z , κ W W , Zu L , Zu R , Zd L , Zd R , and W u L . The ATLAS search for h → W W * reported in Ref. [42] considers the VBF-enriched category in which the detection of two jets consistent with VBF kinematics is required. The expected yields in this category are reported in Table VII of Ref. [42]. After the final selection cuts at 8 TeV with 20.3 fb −1 of integrated luminosity, the expected number of Higgs VBF events in the SM is 4.7 (compared to 5.5 background events) in the eμ sample. Rescaling the number of expected events with the expected HL-LHC luminosity (3000 fb −1 ) and cross section, we expect about 2000 SM Higgs VBF events to be collected by each experiment. In the following, we make a brave approximation and neglect any background events in the fit and assume that the HL-LHC will observe a total of 2000 events compatible with the SM expectations.
As anticipated, a key point to be addressed for a consistent extraction of the PO is the validity of the momentum expansion. In order to control such an expansion, we set an upper cut on the p T of the leading VBF-tagged jet. The momentum expansion of the form factors in Eq. (6) only makes sense if the higher-order terms in q 2 1,2 are suppressed. This requirement leads to the consistency condition, where q 2 max is the largest momentum transfer in the process. A priori we do not know the size of X f or, equivalently, the effective scale of new physics. However, a posteriori we can verify by means of Eq. (20) if we are allowed to truncate the momentum expansion to the first non-trivial terms. In practice, setting a cut-off on p T we implicitly define a value of −q 2 max . Extracting the X f for p T,j < (p T,j ) max ≈ −q 2 max we can check if Eq. (20) is satisfied. Ideally, the experimental collaborations should perform the extraction of the X f for different values of ( p T,j ) max optimizing the range according to the results obtained. In the following exercise we set ( p T,j ) max = 600 GeV which, a posteriori, will turn out to be a good choice in absence of any sizable deviations from the SM.
In our analysis we choose the binning in the double differential distributions in the two VBF-tagged jet p T as {30 − 100 − 200 − 300 − 400 − 600} GeV. We use the UFO implementation of the Higgs PO in the Sherpa Monte Carlo generator [34,41] to simulate VBF Higgs events over the relevant PO parameter space in proton-proton collisions at 13 TeV c.m. energy. Here we employ the VBF selection cuts as listed in Eq. (17) with the additional requirement η j 1 j 2 > 3. We verified that the results of the fit are independent on the precise value of this last cut. Renormalization Analyzing the simulation output, we find expressions for the number of expected events in each bin as a quadratic polynomial in the PO: where a is a label for each bin. Assuming that the HL-LHC "would-be-measured" distribution is SM-like and describing the number of events in each bin with a Poisson distribution, we construct a global likelihood L and evaluate the best-fit point from the maximum of the likelihood. We then define the test statistic, χ 2 = −2 log(L/L max ), as a function of the seven PO. For more details of the statistical analysis see "Appendix".
In Fig. 5, we show in red the 1σ ( χ 2 ≤ 1) and 2σ ( χ 2 ≤ 4) bounds for each PO, while profiling over all the others. The expected uncertainty on the κ Z Z,W W is rather large (with a loosely bounded direction: δκ Z Z ≈ −3δκ W W ), however, in a global fit to all Higgs data, these PO are expected to be much more precisely constrained from h → 4 , 2 2ν decays. The most important conclusion of this analysis is that at the HL-LHC all five production PO can be constrained at the percent level. In the following we test the robustness of this conclusion.
The likelihood obtained from the PO fit is highly non-Gaussian, which is mainly due to the fact that Eq. (21) is quadratic in the PO, and thus the χ 2 is approximately a quartic polynomial. This implies that using the Gaussian approximation to obtain the 1σ uncertainties from an expansion around the minimum overestimates these errors (compare with the 1σ intervals of Fig. 5): In order to assess if these bounds simply come from the information of the total rate, which in a complete analysis depends on the decay parameters and the total Higgs decay width, or it indeed stems from the shape analysis, we introduce a new parameter μ as an overall rescaling of the number of events in all bins, N ev a → μN ev a . We then perform the same fit as above with this extra parameter and subsequently profile over it. 6 As a result, κ Z Z and κ W W become unconstrained but the constraints on the contact terms do not change qualitatively. We thus conclude that their bounds do come from the shape information, i.e. the normalized distributionF( p T j 1 , p T j 2 ).
Furthermore, we have checked that the uncertainties on the entries of the X a matrices, due to the finite statistics of our Monte Carlo simulations, do not impact the fit results. Details of this analysis are reported in "Appendix". The approach sketched there can also be used to estimate the uncertainty of our result caused by missing higher-order theory corrections, most notably NLO electroweak effects. As anticipated, the latter can exceed the 10 % level in VBF [31,32]; however, the largest contributions are due to factorizable corrections (EW Sudakov logarithms and soft QED radiation) that can be re-absorbed by a redefinition of the PO. From the results in Ref. [43] for the related process e + e − → ννh we estimate non-factorizable NLO electroweak corrections to barely reach 10 % in some dedicated corners of the phase space (being typically well below such values in most of the phase space). To be conservative, we assign uncorrelated relative errors of 10 % in each element of the matrices X a , by introducing appropriate nuisance parameters, and redo the fit. Profiling over these nuisance parameters, in the Gaussian approximation, we find the following 1σ uncertainties for the PO: κ Z Z Fig. 6 Allowed deviations in the distribution of the leading-jet p T by varying the PO within the −2 log L/L max < 4 (2σ ) region obtained after the VBF fit. In the left plot we show the absolute number of events in each bin, while in the right one we show the normalized distribution with respect to the total number of events and the bin width = 0.94, κ W W = 0.31, Zu L = 0.022, Zu R = 0.027, Zd L = 0.033, Zd R = 0.055, and W u L = 0.009. Interestingly, comparing these with the Gaussian errors shown above, we conclude that the estimated sensitivity does not worsen significantly, indicating that statistical errors will still dominate. It is worth noting that the theoretical uncertainties are more relevant for the determination of κ Z Z and κ W W and less relevant for the contact terms PO. Now that we have obtained the constraint on the PO, we can a posteriori check the consistency condition of the analysis, namely, that we are in the regime of small deviations from the SM prediction. In Fig. 6, we show the envelope of the allowed deviations in the leading-jet p T distribution, obtained by varying the PO inside the 2σ region. As can be seen, the size of the distribution is well constrained up to 400 GeV. Equivalently, using | X f | 0.01 to check the consistency condition (20), we find 0.01×(600 GeV) 2 /m 2 Z 1, suggesting that we have performed an analysis in a kinematical region where the momentum expansion is indeed reliable.

VH kinematics
Higgs production in association with a W or Z boson are respectively the third and fourth most important Higgsproduction processes in the SM, by total cross section. Combined with VBF studies, they offer complementary handles to limit and disentangle the various Higgs PO. Due the lower cross sections, so far these processes are mainly studied in the highest-rate Higgs decay channels, such as h → bb [44][45][46][47] and h → W W * [48][49][50][51]. The drawback of these channels are large backgrounds, which are overwhelming in the bb case and of the same order as the signal in the W W * channels. In the following we skip over the challenges and the difficulties due to the presence of large backgrounds in these dominant modes, focusing only on V + h decay channels with a good S/B ratio (which should become accessible at the HL-LHC).
In those channels we analyze the prospects for the extraction of the corresponding production PO.
An important improvement for future studies of these channels with the much higher luminosity that will be available, can be obtained scrutinizing differential distributions in specific kinematical variables. In Sect. 2.2 we showed that with this respect the (not always measurable) invariant mass of the V h system is the most important observable in this process, since the form factors directly depend on it. In channels where the invariant mass m V h cannot be reconstructed due to the presence of neutrinos, an accessible kinematical proxy exhibiting a sizable correlation with q 2 is given by the transverse momentum of the vector-boson p T,V or, equivalently, that of the Higgs, as can be seen in the Fig. 7; see also Ref. [52]. Even though this correlation is not as good as the one between the jet p T and the momentum transfer in the VBF Higgs-production channel, a measurement of the vector-boson (or Higgs) p T spectrum would still offer important information on the underlying structure of the form factors appearing in Eq. (9), namely F q i Z L (q 2 ) or G q i j W L (q 2 ); see also Ref. [53]. The invariant mass of the V h system is given by In the c.m. frame, we have p V = (E V , p T , p z ) and p h = (E h , − p T , − p z ) and For p z = 0 this equation gives the minimum q 2 for a given p T , which can be seen as the left edge of the distributions in Fig. 7. This is already a valuable information, especially to address the validity of the momentum expansion. For example the boosted Higgs regime utilized in many bb analyses implies a potentially dangerous lower cut-off on q 2 : here a bin with p T > 300 GeV implies q 2 630 GeV, which might be a problem for the validity of the momentum expansion.
In the W h process, for a leptonic W boson decay, the p T,W can not be reconstructed independently of the Higgs decay channel. It is tempting to consider the p T of the charged  lepton from the W decay as correlated with the W h invariant mass. However, we checked explicitly that any correlation is washed out by the decay.

NLO QCD corrections in VH
At the inclusive and exclusive level QCD corrections to VH processes are well under control [26,27,54]. The dominant QCD corrections of Drell-Yan-like type are known fully differentially up to NNLO [55][56][57] and on the inclusive level amount to about 30 % with respect to the LO predictions for both W h and Zh. Remaining scale uncertainties are at the level of a few percent.
In Fig. 8 we illustrate the NLO QCD corrections to Zh in the SM looking at differential distributions in p T,Z and m Zh , while the qualitative picture is very similar for W h. The employed setup is as detailed already in Sect. 3.2, while here we do not apply any phase-space cuts. Although the natural scale choice for VH clearly is μ 0 = Q = ( p h + p Z ) 2 , here we employ a scale μ 0 = H T /2. With this scale choice the resulting differential distributions (to be utilized in the Higgs PO fit) are almost free of shape effects due to higher-order QCD corrections. A study of a similar stabilization including deformations in the Higgs PO framework will be performed in the near future.
In the case of Zh besides Drell-Yan-like production there are loop-induced contributions in gg → Zh mediated by heavy quark loops, which in particular become important in the boosted regime with p T,H > 200 GeV [58,59]. Besides QCD corrections also EW corrections give relevant contributions and shape effects to VH processes due to Sudakov logarithms at large energies. They are known at NLO EW [60,61] and decrease the LO predictions by about 10 % for p T,Z = 300 GeV and by about 15 % for p T,W = 300 GeV. We stress that, as in the VBF case, the dominant NLO EW effects are factorizable corrections which can be re-absorbed into a redefinition of the PO.

Prospects for the Higgs PO in Zh at the HL-LHC
In order to estimate the reach of the HL-LHC, at 13 TeV and 3000 fb −1 of integrated luminosity, for measuring the Higgs PO in Zh production, we consider the all-leptonic channel Z → 2 , h → 2 2ν. The 8 TeV ATLAS search in this channel [51] estimated 0.43 signal events with 20.3 fb −1 (Table X of [51]). By rescaling the production cross section and the luminosity up to the HL-LHC we estimate approximately ∼130 signal events at the SM rate.
Assuming a sample of this size we perform a fit of the p T distribution of the Z boson. In order to control the validity of the momentum expansion we apply an upper cut of p max T = 280 GeV, which corresponds approximately to q 2 ≈ 600 GeV (see Fig. 7). We bin the p T,Z distribution as where a denotes again the label of each bin. We assume the number of events for each bin to follow a Poisson distribution and we build the likelihood L(κ) as a function of the five PO listed above. The best-fit point is defined by L max and we determine χ 2 = −2 log L/L max . In Fig. 5 we show the resulting 1σ (2σ ) intervals for each PO with solid (dashed) blue lines, when all other PO are profiled. The expected bounds obtained in the Zh channel are comparable in strength with the ones obtained in the VBF channel. In Fig. 9 we illustrate the 2σ allowed deviation of the p T,Z distribution A fit based on a binning of the Zh invariant mass spectrum provides very similar errors as those shown in Fig. 5. Again Gaussian errors obtained by expanding the likelihood as a quadratic function around the minimum overestimates the errors compared to the ones shown in Fig. 5, although here not as badly as in the VBF case: Zh : σ Gauss quad (κ Z Z , Zu L , Zu R , Zd L , Zd R ) = (0.085, 0.012, 0.014, 0.013, 0.019).
By multiplying the number of events in each bin by an overall rate modifier μ, as done above for the VBF analysis, and profiling over this parameter, we find κ Z Z being unconstrained but the 1σ errors on the contact terms, in the Gaussian approximation, are exactly the same as the ones before. This clearly implies that the bounds on the contact terms arise from the shape information, and not from the rate.

Prospects for the Higgs PO in W h at the HL-LHC
In the case of W h production, in all the channels used for the Run-1 analysis, the signal manifests itself as a small excess over a large (dominating) background; see e.g. Ref. [51]. A detailed analysis for such processes should be performed evaluating carefully the backgrounds, which is beyond the scope of this work. However, given the high luminosity we are looking at, the golden channel h → 4 , W → ν becomes an interesting viable possibility. It has been estimated by ATLAS that 67 signal SM events will be present with 3000 fb −1 of integrated luminosity [62]. We have thus decided to analyze the prospects of this clean channel only, to constrain the (κ W W , W u L ) PO, with an analogous likelihood analysis as those performed for the Zh and VBF channels. We have studied in particular the p T,H distribution, as reference observable, applying the same binning and upper cut as in the Zh analysis discussed above. In Fig. 5 we show the resulting 1σ (2σ ) intervals for each PO with solid (dashed) green lines, when the other PO is profiled. In this case the Gaussian approximation works well and provides the following 1σ errors: Upon introducing a total rate modifier μ, as done for the previous channels, the bound on κ W W vanishes when μ is profiled. However, the constraint on the contact term PO W u L remains unchanged, implying that also in this case the bound arises from the shape of the p T distribution.
We conclude the last two phenomenological sections stressing that we have performed simplified estimates of the HL-LHC sensitivity on the contact term PO by separately considering a limited set of collider signatures. It is reasonable to expect that, including all possible signatures and performing a global fit, the sensitivity can significantly improve. However, such a global analysis should also consider the effect of backgrounds, neglected in this study.

Validity of the momentum expansion
The momentum expansion at the basis of the PO decomposition is an expansion on kinematical variables that are experimentally accessible. As such, the radius of convergence of this expansion can be checked, a posteriori, by means of experimental data. In particular, as pointed out in Sect. 3.3, a crucial check is represented by the consistency condition (20), where q 2 max is controlled by ( p T,j ) max in VBF and m V h in VH (or, less efficiently, by p T,Z and p T,H in VH). More generally, the high-momentum behaviors of d 2 σ/d p T,j 1 d p T,j 2 (VBF) and dσ/dm V h (VH) provide a direct probe of the validity of the momentum expansion, or the absence of nearby NP poles.
Besides these direct probes of the high-momentum behavior of the cross sections, a further check to assess the validity of the momentum expansion is obtained comparing the fit performed including the full quadratic dependence of N ev a on the PO, with a fit in which the N ev a are linearized in δκ X ≡ κ X − κ SM X and X . The idea behind this procedure is that the quadratic corrections to physical observable in δκ X and X are formally of the same order as the interference of the first neglected term in Eq. (6) with the leading SM contribution.
If the two fits (linear vs. quadratic) provide similar results, one can safely conclude that the terms neglected in the PO decomposition are indeed subleading. In principle, if the two fits yield significantly different results, the difference might be used to estimate the uncertainty due to the neglected higher-order terms in the momentum expansion. In practice, as will be illustrated below, this estimate turns out to be rather pessimistic and often an overestimate of the uncertainty on the PO.
To access the feasibility of this check, we perform a linear fit for the VBF Higgs production closely following the procedure described in Sect. 3.3. The results obtained in the Gaussian approximation are: Comparing those results with Eq. (22), we conclude that the bounds on the contact terms in the linearized case are significantly weaker (typically one order of magnitude less stringent) than those obtained in the quadratic fit. Similar results are obtained for the Zh analysis, while only in the W h case the two fits give comparable results: Given the events we have simulated are obtained using SMlike distributions, we cannot attribute this large difference from a possible breakdown of the momentum expansion in the underlying distribution. We dedicate the rest of this section to investigate in more detail the origin of the mismatch and how to address it. The most likely explanation for the large difference between linear and quadratic fits reported above is the fact that in the linear fit only a few linear combinations of the PO enter the observables, thus reducing the number of independent constraints one can get. This fact, coupled to the large number of free parameters in VBF and Zh, could explain the loose constraints obtained in the linear fit. If this was true, we should find that in simple models with less parameters the linear and quadratic fit should agree.
To check if the constraints obtained on the contact terms can, in fact, be used to bound explicit new physics scenarios, we employ a simple toy model. To this end, we extend the SM with a new neutral vector boson, Z , coupled to specific fermion currents (to be defined below) and to the Higgs, such that it contributes to VBF and VH (or better Zh) production. Since the goal of this section is to examine the validity of the momentum expansion with an explicit new physics example, we ignore all other phenomenological constraints on such a model (for example, electroweak precision tests, direct searches, etc.). 7 One the one hand, we compute the bounds on the mass and couplings of this new state from the analysis of the double differential p T distribution in VBF Higgs production (and the p Z T distribution in Zh). On the other hand, we integrate out the heavy Z and match to the Higgs PO framework. Finally, we compare the bounds in the full model with the ones obtained from the Higgs PO fit.
To be more specific, we consider a Z which contributes to the form factor F f f Such a contribution could arise, for example, from the following interaction terms: where all the fields are canonically normalized and in the mass basis. Using FeynRules [65] (package version 1.6.16) we obtain an UFO [11] representation of this Z model and perform exactly the same analysis previously applied to the PO for VBF and VH production. This allows us to derive bounds on the combination of couplings g f ≡ g H g f Z for a set of benchmark Z masses, M Z . In this simple model the Z only decays to a pair of fermions as well in Z + h. The corresponding partial decay widths, assuming the Z is much heavier than the daughter particles, are where N c is the number of colors. In order to simplify the analysis, we assume that the Z is a narrow resonance ( Z M Z ). This allows one to interpret bounds from the VBF and VH analyses in terms of the g f parameters. Using the above relations, we have checked that this condition is satisfied for the benchmark scenarios we consider in the following. Expanding the form factor from Eq. (29) for q 2 1 M 2 Z and Z M Z and keeping only the leading deviation from the SM, we find

Effect of the Z in VBF
We consider the case where the Z couples to both the down and the up right-handed quarks, with two independent couplings, g d R Z and g u R Z . In addition, we fix the Z mass to two benchmarks values: (a) 700 GeV and (b) 2000 GeV. The main results of the analysis are shown in Fig. 10.
On the one hand, we perform a fit to the Higgs PO Zu R and Zd R , while fixing all other PO to zero, and translate this bound on the relevant parameter space of the Z model, namely the {g d R , g u R } plane. We report the results of the fit obtained with full quadratic dependence on the PO, as well as the results in which N ev is linearized in δκ X and X . In both cases, 95 % CL bounds are obtained by requiring −2 log L/L max ≤ 5.99. On the other hand, using exactly the same binning and statistical treatment, we directly fit the Z model parameters.
Comparing the two methods we conclude: (i) for both masses the quadratic PO fit provides a reasonable approximation of the model fit, while the linear fit largely overestimates the errors; (ii) the PO fit performs better for M Z = 2000 GeV than for M Z = 700 GeV, as expected from the momentum expansion validity arguments (we recall that we set the cut p T,j < 600 GeV); however, also for M Z = 700 GeV the quadratic fit does provide a fair approximation to the model fit. In particular, in this case we see that the bound from the PO fit is stronger than in the model, which can be understood by the fact that in VBF the Z is exchanged in the t-channel, and therefore its main effect is to reduce the amplitude for high values of q 2 .

Effect of the Z in Zh
In order to assess the validity of the momentum expansion in associated production, it is convenient to look first at the underlying partonic cross section. In Fig. 11 we show the partonic cross section dd → Zh, as a function of the Zh invariant mass, for the two benchmark points of Z model introduced above.
Both benchmark points have been chosen such that they generate the same contact term when the Z is integrated out, Zd R = 1.68 × 10 −2 , which is within the 2σ bound of our PO fit. The width of the Z has been fixed to 100 GeV and 200 GeV for the light and heavy scenario, respectively. Using Eq. (31) and assuming no other decay mode is present, this corresponds to g H 0.097 (3.0) in the light (heavy) scenario. We have checked that our conclusions do no change by varying the total width, as long as the condition Z M Z is satisfied. As expected, in the light scenario the cross section in the full model strongly deviates from the PO one well before the 600 GeV cut-off imposed in the fit, implying that our PO fit is not reliable in this case. On the other hand, the scenario with a heavy and strongly coupled Z shows a very good agreement with the full PO analysis up to ∼1 TeV, i.e. well above the UV cut-off of our analysis, implying that the analysis can be safely applied to such scenarios, and that it could be even improved by setting a slightly higher cut-off. In both cases, from Fig. 11 is clear that the linearized dependence on the PO is not sufficient to describe the cross section, even for energies much smaller than the Z mass.
From this analysis we can anticipate the results of a comparison of various fits of Zh data, i.e. full model fit vs. PO fits using quadratic and linear dependence, as already done in the VBF case. In Fig. 12 we show the results of such fits. We stress that in all cases the analysis was exactly the same: we have analyzed the p Z T distribution up to 280 GeV, employing always the same binning (as discussed in Sect. 4.3). The solid red line represents the 95 % CL bound in the full model while the solid (dashed) blue line shows the bound obtained from the PO fit with quadratic (linear) dependence.
The distributions in Fig. 11 allow a straightforward interpretation of these results. In the heavy-Z case, the full quadratic expansion in the Higgs PO describes very well the m Zh distribution before the cut-off of 600 GeV, while keeping only the linear dependence underestimates the new physics contribution. It is thus expected that in this case the bound will be much worse. In the light-Z case, both From this illustrative toy-model example we can draw the following general conclusion with respect to the validity of the PO expansion: for underlying models that respect the momentum expansion, hence for models where the PO extracted from data satisfy, a posteriori, the consistency condition (20), the quadratic fit provides more reliable and thus more useful constraint on the PO. In such models the difference between quadratic and linear fit represents a large overestimate of the errors.
However, the situation is more involved for models with low-scale new physics. The latter should manifest by anomalously large values of the PO, or sizable differences in the fits performed with different upper p T cuts. In such cases the quadratic fit is likely to provide a useful constraint, especially for the class of models with a strong correlation between linear and quadratic terms in the momentum expansion (as the simple Z model discussed above). Still, for low-scale new physics we cannot exclude more complicated scenarios where new model parameters appearing at higher order in the momentum expansion wash-out an apparent small error on the PO from the quadratic fit. In such cases only the results of the linear fit (with a properly low p T cut) would provide an unbiased constraint on the model.
In view of these arguments, we encourage the experimental collaborations to report the results of both linear and quadratic fits, as well as to perform such fits using different p T cuts. Table 1 Summary of the "production PO", namely the PO appearing in VBF and VH in addition to those already present in Higgs decays (classified in Ref. [1]). In the second column we show the independent PO needed for a given set of amplitudes, assuming both CP invariance and U (2) 3 flavor symmetry. The additional variables needed if we relax these symmetry hypotheses are reported in the third and fourth columns. In the bottom row we show the independent PO needed for a combined description of VBF and VH under the hypothesis of custodial symmetry. The number of independent PO range from 12 (sum of the first two lines) to 4 (bottom row, second column) Amplitudes/processes U (2) 3

Conclusions
Higgs physics is entering the era of precision measurements: future high-statistics data will allow us not only to determine the overall signal strengths of production and decay processes relative to the SM, but also to perform detailed kinematical studies. In this perspective, an accurate and sufficiently general parameterization of possible NP effects in such distributions is needed. In this paper we have shown how this goal can be achieved in the case of VBF and VH production, generalizing the concept of Higgs PO already introduced in Higgs decays.
As summarized in Table 1, the number of additional PO appearing in all VBF and VH production amplitudes is manageable. In particular, assuming CP invariance, flavor and custodial symmetry, only four new PO should be added to the set of seven PO appearing in h → 4 , 2 2ν, 2 γ , 2γ in the same symmetry limit [1]. This opens the possibility of precise global determinations of the PO from combined analyses of production and decay modes, already starting from the next LHC runs.
As extensively illustrated in Sects. 3 and 4, the key aspects of VBF and VH is the possibility of exploring sizable momentum transfers in the Green functions of Eq. (1). On the one hand, this maximizes the sensitivity of such processes to PO that are hardly accessible in Higgs decays. On the other hand, it allows us to test the momentum expansion that is intrinsic in the PO decomposition as well as in any EFT approach to physics beyond the SM. Key ingredients to reach both of these goals are precise differential measurements of d 2 σ/d p T,j 1 d p T,j 2 in VBF and dσ/dm V h in VH (or appropriate proxies such as p T,H and p T,Z ). We thus encourage the experimental collaborations to directly report such differential distributions, especially in the kinematical regions corresponding to high momentum transfer.
As far as the PO fits in VBF and VH are concerned, we suggest to perform them setting a maximal cut on p T,j and m V h , to ensure (and verify a posteriori) the validity of the momentum expansion. As illustrated by matching the PO framework to simplified dynamical NP models, it is also important to report the results of fits using both linearized and quadratic expressions for the cross sections in terms of PO. According to our preliminary estimates, the production PO could be measured at the percent level at the HL-LHC (in the case of maximal flavor symmetry, without the need of imposing custodial symmetry). This level would be sufficient to constrain (or find evidence of) a wide class of explicit NP models and, among other things, to perform non-trivial tests of the relations between electroweak observables and Higgs PO expected in the SMEFT.

Appendix: Details of the statistical analysis
In this appendix we provide details of the statistical analysis used to derive the projected sensitivity on the PO. The first step of such an analysis is to compute, by means of Monte Carlo simulation, the signal yield in each bin as a quadratic polynomial in the PO: where a labels a given bin. For example, the total cross section for VBF Higgs production at 13 TeV, applying the cuts defined in Sect. 3.3, is given by After deriving similar expression for each bin, we perform a profile likelihood fit to a binned histogram distribution. The likelihood function takes the form where N exp a denotes the number of projected-observed events in a bin with label a (which we take to be SM-like), X a i j are the coefficients of the X a matrix as obtained from our simulation and X a i j are the uncertainties associated with these coefficients. These uncertainties are determined from a Poisson distribution in the number of events in each bin and from a normal distributions in the nuisance parameters X a i j . We first minimize the above function with respect to the PO (κ) and nuisance parameters (X a i j ) and then expand the function around the best fit point up to second order where dots represent terms that involve the nuisance parameters as well. Here, V i j = σ i ρ i j σ j where σ i and ρ i j are the uncertainties and correlation coefficients, respectively. We refer to this method as the Gaussian approximation, and use it to study the impact due to Monte Carlo uncertainties, as well as due to missing higher-order corrections on the fit results. On the other hand, the results shown in Fig. 5 are obtained by setting the error on the nuisance parameters to zero, and for each PO, profiling over all the others.