A sensitivity study of VBS and diboson WW to dimension-6 EFT operators at the LHC

We present a parton-level study of electro-weak production of vector-boson pairs at the Large Hadron Collider, establishing the sensitivity to a set of dimension-six operators in the Standard Model Effective Field Theory (SMEFT). Different final states are statistically combined, and we discuss how the orthogonality and interdependence of different analyses must be considered to obtain the most stringent constraints. The main novelties of our study are the inclusion of SMEFT effects in non-resonant diagrams and in irreducible QCD backgrounds, and an exhaustive template analysis of optimal observables for each operator and process considered. We also assess for the first time the sensitivity of vector-boson-scattering searches in semileptonic final states.


Introduction
While the CERN Large Hadron Collider (LHC) already collected an unprecedented amount of data, no significant deviations from the Standard Model (SM) predictions have been observed so far. In the upcoming Runs, the LHC experiments will be able to look for physics Beyond the SM (BSM) in new ways, complementing the quest for patent signatures, such as new resonances, with searches for small deviations with respect to the SM expectations. This will be made possible by the increasing statistical precision granted by the available data sets, together with the improved understanding of the detector behaviour and the application of advanced data-analysis techniques based on state-of-the-art machine learning algorithms.
The primary theoretical tool for BSM searches in precision measurements is the Standard Model Effective Field Theory (SMEFT) [1][2][3] (see [4] for a review). The SMEFT has been developed extensively in the past decade. State-of-the-art analyses are based on Higgs, diboson and top quark measurements and constrain about 20-30 SMEFT parameters simultaneously [5][6][7][8][9][10][11][12]. Going forward, the scope of these studies will be extended in several directions. For instance, there is interest in relaxing CP and flavour assumptions [13], as well as in incorporating information from B-meson observables [14][15][16]. The quality of the fit ingredients is also expected to improve in the future: on the theoretical side, more accurate SMEFT predictions are becoming available, that include 1-loop QCD or electro-weak (EW) corrections. On the experimental side, current measurements will improve substantially, achieving higher precision and allowing the extraction of differential information with higher resolution. Several processes will also become accessible at the LHC for the first time.
Vector Boson Scattering (VBS) ones are among the most interesting signatures in this perspective [17]. The analysis of Run-II data by the ATLAS and CMS collaborations has recently led to the observation of VBS in the same-sign and opposite-sign WW [18][19][20], WZ [21,22] and ZZ [23,24] final states, and to strong evidence in the semileptonic WV final state [25]. VBS processes with a photon, a heavy vector boson and two jets were recently observed as well [26][27][28]. This broad class of processes is mainly statistically limited at the LHC. Therefore the increased luminosity of its upcoming Runs will benefit their measurements in a particularly significant way, allowing the extraction of new constraints on the EW sector, that are complementary to those from Higgs and diboson measurements. Specifically, VBS processes provide tree-level sensitivity to effective operators inducing modifications of triple (TGC) and quartic gauge couplings (QGC), as well as of Higgs-gauge couplings away from the Higgs mass-shell, and even to contact interactions among four light quarks.
Mainly motivated by the QGC sensitivity, previous studies of EFT effects in VBS were often restricted to dimension-8 operators that generate genuine QGC corrections while leaving TGCs unaffected [29], see e.g. Refs. [30][31][32][33]. The impact of dimension-6 operators was explored systematically only recently, driven by the interest in incorporating these measurements into global SMEFT analyses. Even though VBS measurements typically yield weaker constraints -in absolute terms -compared to EW precision observables, Higgs or inclusive diboson measurements, they can play a significant role in global analyses by constraining new directions in the parameter space and providing a link between the EW, Higgs and four-quarks sectors. Dedicated studies at dimension 6 were carried out for the ZZ [34] (see also [35]) and same-sign WW [36] channels. A first global analysis, that combines VBS and diboson constraints in several channels, was presented recently in [10], demonstrating that, while at present VBS has a visible but small impact in global EFT fits, in the future it will play a more significant role.
In this work we examine the individual and combined sensitivity of five different VBS channels to a set of 14 dimension-6 operators in the Warsaw basis [2]. The primary goals of our study are to explore the complementarity among different VBS channels and to determine their relative projected sensitivities to different classes of effective operators. We also aim at assessing the relevance of SMEFT contributions that are quadratic in the Wilson coefficients, and of SMEFT corrections to the irreducible QCD backgrounds.
We treat VBS processes as 2 → 6 scatterings, retaining EFT contributions to non-resonant diagrams. As these processes exhibit a great kinematic complexity, we perform a comparative study of several kinematic variables for each channel, to determine which observables are most sensitive to each dimension-6 operator. The expected performance of the VBS channels is compared to that of diboson processes, which are traditionally the main playground for EFT analyses at dimension-6, due to their large production cross-sections and reasonable signalover-background ratios. This is achieved by including inclusive WW diboson production [37,38] in our analysis.
For simplicity, we work at parton-level and at leading-order (LO) in QCD, and we neglect all backgrounds other than diboson-plus-jets. 1 Due to these simplifications, that prevent us from reproducing a fully realistic analysis, we refrain from confronting our results with current measurements. Instead, we extract expected limits on the Wilson coefficients for integrated luminosities of 100 fb −1 or larger, that should be interpreted as indicative of the maximal sensitivity of these processes to EFT effects.
The manuscript is organised as follows: the theory framework is defined in Section 2. The adopted strategy for the numerical and statistical analysis is presented in Sec. 3.2, while the specific properties of the individual processes are discussed in Sec. 4. We present our results in Sec. 5 and in Sec. 6 we conclude.

Effective Field Theory framework
We consider the SMEFT truncated at dimension-six: where Λ is the cutoff scale of the EFT, c α are the Wilson coefficients, and the index α runs over the labels of a complete and non-redundant set of operators that are invariant under the SM gauge symmetries. Here we work in the so-called Warsaw basis [2] and the Lagrangian is defined as in the SMEFTsim package implementation [39,40]. In particular, we require the fermionic operators to be invariant (up to insertions of the Yukawa couplings) under a U (3) 5 flavour symmetry. We also assume CP conservation and take the Fermi constant (G F ) and the W and Z boson pole masses (m W , m Z ) as input quantities for the EW sector [41,42]. The SU (2) Higgs doublet is denoted by H, the gauge field strengths associated to the SU (2) and U (1) symmetries by W i µν and B µν respectively, the left-handed lepton and quark 1 In most cases, diboson-plus-jets constitutes the dominant background. The main exception is ZV jj production in the semileptonic channel, for which the dominant background is Z+jets. Further (reducible) background sources can come from particle mis-identification or similar detector effects, and cannot be reproduced in a parton-level analysis.
The subset of Warsaw basis operators considered in this work. Repeated indices are understood to be summed over. p, r are flavour indices, and a U (3) 5 -invariant flavour structure is assumed.
doublets by l, q and the right-handed quark and charged-lepton fields by u, d, e. The SU (2) indices are indicated with i, j, k and the Pauli matrices by σ i . Flavour indices are indicated with p, r. For further notational conventions we refer the reader to Ref. [40].
For our numerical analysis, we consider only the subset of 14 operators given in Table 1. It contains all the operators that enter via modifications of the EW input quantities (Q Hl , Q (1) ll , Q HD , Q HW B ), plus a set of dimension-6 operators that give significant contributions to all VBS processes, once the experimental selection cuts are applied. These are mainly induced via modifications of V f f (Q  Table 1 represents a convenient selection of operators for the purposes of this work, as it allows to examine all the mentioned categories of SMEFT effects and to explore the complementarity between VBS processes in constraining EFT parameters, while at the same time avoiding a very high-dimensional fit space. Fermionic operators with right-handed fermions and bosonic operators such as Q HB , that only enter a subset of VBS processes, as well as contact interactions between two quarks and two leptons, would need to be added for an exhaustive global analysis. These would amount to about 20 extra degrees of freedom, most of which are not expected to introduce new significant features to the fit. In this sense, we consider the set in Table 1 adequate for a study of the sensitivity of VBS processes to EFT effects and we leave a more complete analysis for future work.
Working at order Λ −2 , a generic scattering amplitude has the form: where A SM is the SM amplitude and A Qα is the total amplitude obtained with one insertion of the operator Q α . The latter scales linearly with c α /Λ 2 , and this dependence has been made explicit in Eq. (2.2).
As a consequence, the expected number of events in a given phase-space region scales with the Wilson coefficients as: where in Eq. (2.3) both α and β run over all relevant indices, and for α = β the last contribution reduces to (c 2 α /Λ 4 )|A Qα | 2 . The final result is therefore the sum of the SM prediction, a term that is linear in the Wilson coefficients and stems from the interference between SM and BSM amplitudes, and a pure BSM contribution, that is quadratic in the Wilson coefficients and is usually referred to as "quadratic" term. In Eq. (2.4) the latter has been split into individual and mixed quadratic contributions for later convenience.
The result in Eq. (2.4) applies both to integrated observables and bin-by-bin to differential ones. For a given observable, the quantities N SM , N int α , N quad α , N mix α,β can be estimated numerically. For n operators contributing, there are n independent linear terms, n individual squares and n(n − 1)/2 mixed ones. Therefore, the whole set of EFT contributions can be determined by computing N at a total of n(n + 3)/2 (119 for n = 14) independent points in the parameter space {c α }. In this work, these estimates are performed via Monte Carlo event simulations, as described in the next section.
Note that insertions of a SMEFT operator into the scattering amplitude can take place both in vertices and in the propagators of unstable particles, via corrections to their masses and decay widths. As all masses are taken to be input quantities, only corrections to the decay widths of the W, Z and Higgs bosons are relevant to our calculation. Width corrections are polynomial in the Wilson coefficients and enter the denominator of the scattering amplitude. In order to recover the form in Eq. (2.2) they are typically Taylor-expanded in c α . The linear corrections obtained in this way are well-defined [43] and can be estimated directly with available Monte Carlo tools [40]. In practice, one has where N int α,vert. is the contribution from vertex insertions while the N int α,δΓ terms come from insertions in the W, Z and H propagators respectively. In the Monte Carlo simulation, each of the terms in (2.5) is estimated separately.
The propagator contributions to N quad α , N mix α,β , on the other hand, are ambiguous. The source of ambiguity is that the quadratic terms defined above provide partial predictions at order Λ −4 , as other contributions of the same order, such as insertions of two EFT operators in the same amplitude, are neglected. While for vertex corrections the distinction between double insertions and "standard quadratics" is well-defined, for propagator corrections it is less clear. For instance, the result obtained expanding first A to linear order in c α and then taking a narrow-width approximation (NWA) on |A Qα | 2 is different from the one obtained taking first the NWA and then expanding the SMEFT corrections to Γ in the (σ prod. Γ part. /Γ) expression. Agreement between the two procedures is recovered if, in the first case, the square of the linear correction to the width (δΓ lin ) 2 is retained in the expanded A Qα .
As, with currently available tools, it is not possible to produce complete simulations to quadratic order in the dimension-6 Wilson coefficients, we choose to omit all propagator corrections from our baseline analyses. Although admittedly not ideal, this solution ensures a consistent comparison between linear and quadratic results and avoids introducing arbitrariness in our predictions. Nevertheless, we do estimate the impact of propagator corrections on the purely linear analysis, as discussed in Sec. 5.1.

Analysis procedure
We consider five distinct VBS channels, characterized by the presence of two well-separated hadronic jets in the final state, plus the inclusive production of opposite-sign W pairs, that serves as a representative of diboson processes. As indicated in Table 2, all VBS channels are defined as 2 → 6 scatterings, while inclusive WW diboson production is defined as a 2 → 4 process. The channels considered are: • W ± W ± +2j: same-sign WW (SSWW) production with two same-sign leptons, missing transverse energy and two jets in the final state; • W + W − +2j: opposite-sign WW (OSWW) production with two opposite-sign leptons, missing transverse energy and two jets in the final state; • W ± Z+2j: WZ production with three charged leptons, missing transverse energy and two jets in the final state; • ZZ+2j: ZZ production with four charged leptons and two jets in the final state; • ZV +2j: ZV production (combining W ± Z and ZZ) with a semileptonic final state, where the Z boson decays into two charged leptons and the V one decays into two quarks, plus two tagging jets; • W + W − : inclusive WW production with two opposite-sign leptons and missing transverse energy. As we work at LO in QCD, this channel has no additional jets.
All Feynman diagrams leading to the relevant final states are retained in our signal definitions, including those with non-resonant topologies. In the fully-leptonic channels, the final states include electrons and muons and, for simplicity, the two vector bosons are required to decay into lepton (or lepton-neutrino) pairs of different flavor. For all VBS channels except SSWW+2j, irreducible QCD backgrounds are present in the SM, that also receive SMEFT corrections. These contributions are taken into account in our study and their role in the SMEFT fit is discussed in Sec. 5 ZZ+2j EW generate p p > e+ e-mu+ mu-j j QCD=0 SMHLOOP=0 ZV+2j EW generate p p > z w+(w-,z) j j QCD=0 SMHLOOP=0, z > l+ l-, w+(w-,z) > j j OSWW+2j QCD generate p p > e+ ve mu-vm˜j j QCD==2 SMHLOOP=0 WZ+2j QCD generate p p > e+ e-mu+ vm j j QCD==2 SMHLOOP=0 ZZ+2j QCD generate p p > e+ e-mu+ mu-j j QCD==2 SMHLOOP=0

Event generation
For each channel we simulate events at parton-and tree-level with Madgraph5_aMC@NLO (v. 2.6.5) [44], interfaced to the SMEFTsim package (v. 3) [39,40]. The syntax used to generate each process is reported in Table 2 for the SM case 2 , together with a specification of how the syntax is modified for the extraction of the SMEFT contributions. The simulations are produced for pp collisions at a centre-of-mass energy of 13 TeV, using the NNLO parton distribution functions provided by the NNPDF collaboration [61], with α S = 0.118 and in the four-flavour scheme (LHAPDF identification code 325500). The renormalization and factorization scale choices are determined by the Monte Carlo generator as the transverse mass of the 2 → 2 scattering after k T clustering. For each process, the SM, interference and quadratic EFT components (Eq. (2.4)) are extracted using two alternative techniques: a) performing an event simulation for each contribution. The SM signal (N SM ), the pure interference (N int α,vert. and N int α,δΓ ) and quadratic term (N quad α ) for individual operators, as well as the mixed quadratic contributions N mix α,β can all be generated directly exploiting the interaction-order syntax in Madgraph5_aMC@NLO and SMEFTsim. b) exploiting the reweighting method in Madgraph5_aMC@NLO [62]. In this case the events are generated only once, for one point in parameter space, and are subsequently reweighted to match a different set of values of the Wilson coefficients. Because this procedure is based on the calculation of the matrix element at fixed phase-space points, the dependence of the weights on the Wilson coefficients is computed exactly.
While N int α , N quad α can be determined directly, isolating the quantities N mix α,β often requires algebraic combinations of results obtained at 2 or 3 different points. Besides being faster, the reweighting technique has the advantage of reducing significantly the statistical uncertainties in this operation. On the other hand, it can give unreliable results in regions of phase space that are poorly populated in the original event generation. For this reason, both techniques were used for estimating and cross-checking the signal dependence on the Wilson coefficients.
Linear corrections stemming from SMEFT insertions in vertices and in the propagators of the W, Z and Higgs bosons were simulated separately, with the syntax indicated in Table 2. Note that propagator corrections have identical shapes for all operators contributing. For instance, for any observable: and similarly for the other N int α,δΓ terms. In this way, all N int α,δΓ can be estimated with one simulation per process and heavy boson (W/Z/H). For later convenience, we report here the numerical expressions of the linearized width corrections (retaining only the operators in Tab. 1): Hq − 1.37c Hl − 0.07c HD + 0.46c HW B , Hl + 1.40c Table 3 summarizes which dimension-six operators contribute to each of the processes under investigation, when they are considered individually. As mentioned above, all operators considered here contribute to all VBS channels, both at interference and quadratic level. However, some of them do not enter certain QCD-induced backgrounds or inclusive diboson (WW) production. Also, certain operators only contribute to non-resonant diagrams, i.e. where the lepton pairs in final state do not come from an intermediate gauge boson. These are indicated with brackets in the table. A priori, contributions from the W, Z and Higgs propagators can introduce a dependence on operators that do not enter via vertices. However, we have verified that this is never the case for the processes and operators considered.

Event selection and analysis strategy
For each process under study, the simulated events are first filtered applying kinematic cuts that reproduce typical acceptance regions for lepton and jet reconstruction in LHC experiments, and isolate a phase space region where VBS production is strongly enhanced over the background. For this study, we refrain from applying any unitarisation procedure or clipping of the high-energy distribution bins. Following a procedure similar to Refs. [63][64][65], we perform a template analysis based on the distributions of the SM, linear and quadratic signal components as a function of a set of final-state-dependent leptonic and partonic observables.
The full list of variables considered is reported in Table 4, along with the selection cuts applied to each final state, and they are defined as follows: • p T,l i , the transverse momentum of the i-th charged lepton (sorted by their p T , from largest to smallest) with respect to the beam axis, η l i its pseudo-rapidity; • m ll , the invariant mass of the two leading-p T leptons (independently of their charge and flavor); • MET, the missing transverse energy, calculated as the transverse component of the vector sum of all neutrino momenta; • p T,j i , the transverse momentum of the i-th outgoing parton (sorted by p T ), η j i , φ j i its pseudo-rapidity and azimuthal angle respectively; • m jj , the invariant mass of the two highest-p T outgoing partons, ∆η jj their pseudorapidity separation and ∆φ jj their angular separation in the azimuthal plane; • ∆R(l, j), the lepton-parton separation in the distance parameter R, defined as the sum in quadrature of φ and η distances.
For the remaining observables we use an analogous notation, so their definitions can be easily inferred. The angular variables Φ planes , θ lW , θ lZ and θ * for the WZ+2j channel are defined in Sec. 4.4 below. For the semileptonic ZV+2j channel, two dijet pairs are formed based on the highest invariant mass as described in Sec. 4.6 and the observables associated to such systems are labelled as "mjj max" or "mjj nomax" with clear connection to the result of the tagging procedure. For each process and EFT operator, only the variable that gives the strongest constraint on the Wilson coefficient (determined a posteriori) is employed in the construction of the likelihood used for the results extraction, as described below.

Likelihood construction
The likelihood function L for the data to match the EFT model is built using Eq. (2.4) to determine the mean of the Poisson statistics (N k (c)) describing the expected number of events surviving the data analysis selections, either in total or in the k-th distribution bin: where c is the vector of free Wilson coefficients in the fit, N k the expected number of events defined as in Eq. (2.4) and n k ≡ N k (0) is the expected number of events in the SM. Where relevant, the QCD SM and EFT components are accounted for in the expression for N k and summed to their EW counterparts. The number of expected events is always normalised to an integrated luminosity of 100 fb −1 , unless otherwise stated. No systematic uncertainty is taken into account, except for a correlated 2% variation between all yields, samples, and  Table 4. Summary table for processes, variables and selections considered in this work. The second column lists the observables examined, while the third one summarises the parton level phase space definition used in this analysis. The last column reports the expected SM event yields of the EW and QCD-induced processes after the analysis selections, for an integrated luminosity of 100 fb −1 . The charged lepton shortcut l stands for electron or muon.
bins, as a proxy of a typical LHC luminosity uncertainty. While actual data analyses will be more severely affected by systematics and theory uncertainties, we refrain from providing estimates for the latter, as they can vary significantly depending on the channel and on the analysis details. This choice is meant to avoid introducing this arbitrariness into the sensitivity estimates, and implies that our results should be interpreted as optimistic estimates of the sensitivity of VBS processes to EFT effects.
For each process considered, only one of the kinematic distributions listed in Table 4 is employed at each time to construct the likelihood. For one-dimensional fits, we choose the variable that gives the strongest limit at 68% confidence level (c.l.) on the Wilson coefficient of interest. Note that, due to this procedure, constraints on different Wilson coefficients are generally derived with different, optimized likelihoods. The ranking procedure is applied separately for individual fits where the quadratic terms in the signal parameterization are retained or discarded. It is also repeated for each two-dimensional fit, retaining the observable that gives the smallest area inside the 68% c.l. contour in the bi-dimensional plane. For the 14dimensional fit, the profiled constraint on each Wilson coefficient is derived from a likelihood that implements, for each process, the same observable as in the corresponding 1D case. The full lists of optimal variables employed for each fit are reported in Tables 5-10 in Appendix A.
The expected sensitivity to the Wilson coefficients is estimated based on the likelihood profile: for single Wilson coefficient scans, the 68% and 95% c.l. intervals in the coefficient estimates are determined by requiring −2∆ log L < 1 and −2∆ log L < 3.84 respectively, being ∆ log L the variation of the likelihood logarithm with respect to its maximum value. For bi-dimensional scans, the intervals are instead −2∆ log L < 2.30, −2∆ log L < 5.99 [66].

Processes considered
In this section we discuss the main properties of the processes considered and their role in the SMEFT analysis. The number of simulated SM events passing the selections, reported in Table 4, provides an indication of the signal cross sections and of the typical signal-tobackground ratios. An example of the impact of dimension-six operators on the kinematic distributions is shown for the SSWW+2j channel (see Sec. 4.2) in Fig. 1. Further representative kinematic distributions are provided for all channels in Appendix B.

Inclusive W + W −
The inclusive W + W − production is well-studied experimentally by the ATLAS and CMS Collaborations, see e.g. [37,38] for the most recent results. The large cross-section of this process makes it an ideal environment for the study of new phenomena at high energies. Indeed, integrated and differential measurements of the WW production are included in most global analyses of dimension-six EFT operators, see Refs. [5,6,[9][10][11][12] for recent examples. This process is also a crucial component of the first combined EFT interpretations of EW/Higgs measurements that were recently published by the ATLAS Collaboration [67,68]   way for larger global EFT fits within the LHC experiments. In our analysis we neglect all backgrounds to this process as, in the phase space defined in Table 4, their yield is approximately a half compared to that of the signal. EFT corrections to the inclusive WW production at the LHC are very well-studied in the literature and NLO QCD corrections are known, see e.g. Refs. [69][70][71][72][73][74]. Here we simulate this process at LO in QCD, where WW only occurs with a qq initial state, via the annihilation of the qq pair or by a t-channel exchange of a parton. Within this approximation the main EFT effects are modifications of triple gauge couplings (c W ) and of couplings of the weak bosons to fermions (c Hl ), as well as input shift corrections (c HD , c HW B , c (1) ll ). For most of these coefficients, the WW production gives more stringent constraints than VBS (in fully-leptonic final states), mainly due to the larger cross section. The only exception is the custodialviolating c HD , that is best probed in processes where the dominant Feynman diagrams involve Z or photon couplings to fermions. c (1) Hl is poorly constrained in this final state, because it only enters non resonant diagrams with a Z boson in s-channel. Higgs (c H , c HW ) and four-fermion operators only enter at 1-loop in QCD, and are therefore neglected here. Distributions for interesting variables for some operators this process is sensitive to, after the selections, are shown in Figure 11.
The analysis of the EW production of a same-sign W boson pair (SSWW+2j) selects the leptonic final state in association with two jets, where the W bosons decay in an electronmuon pair plus the corresponding neutrinos. This process has been observed by the LHC Collaborations in the fully-leptonic final state [18,19,75]. With two leptons of the same charge, moderate MET, and two jets with a large rapidity separation and a large dijet mass, the EW W ± W ± process presents a clean signature in the detector. The LO QCD-induced background, with diagrams with exactly two QCD vertices, is small compared to the EWinduced production and can be kinematically separated from the signal. Thus it can be safely neglected in our analysis. The main source of reducible background for this channel are jetinduced fake leptons, mostly stemming from b quarks in tt events. They can be neglected as their impact is very suppressed in the high-energy tails, where we have most of the EFT sensitivity.
As shown in Fig. 2, this analysis is particularly sensitive to 4-quark operators, that induce very marked shape distortions both in energy and angular variables. At quadratic level, SSWW+2j has the highest discriminating power among their Wilson coefficients, see e.g. Fig. 8. On the other hand, it yields the weakest bounds on c HD , as Z boson contributions are subdominant, and on c HW B : the latter dominantly enters via corrections to the weak mixing angle, to which SSWW has limited sensitivity.
Distributions for interesting variables for some operators this process is sensitive to, after applying the selection cuts, are shown in Figure 12.
Despite possessing the largest cross section among the VBS processes, the purely EW production of a pair of opposite-sign W bosons (OSWW+2j) has been observed experimentally only very recently by the CMS Collaboration [20]. The main challenge in accessing this process is the very large irreducible QCD background (cfr. Tab. 4). We analyze the fully leptonic final state, characterized by the presence of two jets, with high energy and a large separation in pseudo-rapidity, two opposite-charge leptons and missing transverse energy. The same-flavour final states (e + e − νν, µ + µ − νν) are overwhelmed by Drell-Yan events, so here we restrict to different-flavour combinations (e ± µ ∓ νν), that have smaller backgrounds and hence a higher sensitivity. The QCD induced production of W + W − +2j (order α 4 EW α 2 QCD ) is sizeable and taken into account in the analysis. The selections in Table 4 reduce the background contamination, which mainly arises from tt production.
The Vector Boson Fusion Higgs contribution is a subdominant component of the VBS signal, once the selections are applied. Nevertheless, some sensitivity to Higgs operators is retained. In fact, thanks to the presence of s-channel diagrams and to the relatively large cross section, OSWW is the most sensitive channel to anomalous Higgs couplings (c H , c HW ), among those considered in this work. Distributions for interesting variables for some operators this process is sensitive to, after the selections, are shown in Figure 13.

W ± Z+2j
We consider the EW WZ production in the fully leptonic channel plus two jets, with the Z boson decaying into electrons and the W boson decaying into a (µ, ν µ ) pair. One of the advantages of this decay mode is the high purity of the multi-leptonic final state. Additionally, the choice of different flavour decays for the Z and the W bosons provides an efficient discrimination between the two bosons. The WZ+2j channel has a small cross-section but was observed at Run II by both ATLAS and CMS [21,22]. The dominant QCD background for the WZ+2j state is the QCD radiation of partons from an incoming quark or gluon, and it is accounted for in the analysis.
The presence of a single neutrino in the final state allows to kinematically reconstruct its momentum component along the z-axis, by imposing a m W constraint on the invariant mass of the W boson decay products. Once the four-momenta of all the final state particles are known, one can extract the invariant mass of the WZ system m W Z , the angular separations between the vector bosons δη W Z , δφ W Z and the separation between their decay planes Φ planes . Furthermore, knowing the collision centre-of-mass boost along the beam direction, more sophisticated angular observables can be constructed, such as the emission polar angles of the leptons with respect to the direction of the decaying bosons in the rest frame of the latter, θ lW and θ lZ , and the vector bosons emission angle in the centre-of-mass reference frame, θ * .
Among the channels considered, WZ+2j is the one with highest sensitivity to c Hl (together with ZZ+2j). The sensitivity to c HD and c HW B is also quite good: as mentioned above, the explicit presence of a Z boson (and also of photons) in the dominant diagrams enhances the sensitivity to these two operators, compensating for the smaller cross section compared to other VBS channels and even compared to the inclusive WW.
Distributions for interesting variables for some operators this process is sensitive to, after the selections, are shown in Figure 14.

ZZ+2j
The observation of the EW ZZ production was achieved by the ATLAS Collaboration using the full Run-II data set [23], while the CMS Collaboration has recently reported strong evidence for this process [24,76]. This VBS channel has the smallest cross-section among all the processes considered, and is one of the rarest SM processes observed to date. Since the QCD production associated to two jets is a dominant background with respect to the VBS signal, the search for this process is very challenging, despite a very clean experimental signature.
Both the ATLAS and CMS Collaborations use multivariate analyses to isolate the EW signal over the large QCD background, after some fiducial cuts. Here we perform a much simpler analysis using the selections listed in Table 4 and considering exclusively the 2e2µ final state. The variables of interest include, besides those in common with other VBS channels, the invariant mass and total transverse momentum of the lepton pair with the largest transverse momentum m l 1 l 2 , p T,l 1 l 2 , the invariant mass of the four-lepton system m 4l , the total transverse momentum of the same-sign lepton pair p T,e ± µ ± and the transverse momentum of the dilepton system e + e − or µ + µ − with invariant mass closest to m Z , taken as a proxy for p T,Z .
A detailed EFT study of this particular process (restricting to resonant diagrams, which in this case is a good approximation) was presented in Ref. [34]. The sensitivity of ZZ+2j to EFT effects is quite limited due to the low total cross section and small signal-to-background ratio. Nevertheless, this channel is relevant for constraining operators that affect specifically Z interactions, such as c (1) Hl . It is also competitive with the other processes for c HD , c HW B . Distributions for interesting variables for some operators this process is sensitive to, after the selections, are shown in Figure 15.

ZV +2j
Even though the LHC collaborations have mainly studied VBS processes in the fully leptonic final states, semileptonic ZV production also ranks among the processes of greatest interest, as it benefits from the larger branching ratio of V → hadrons (about 67% for W bosons and 70% Z bosons). For instance, ZV+2j was shown to be among the most sensitive channels to dimension-8 EFT effects [77,78].
We consider final states with four jets and a pair of electrons or muons, addressing the challenge to identify properly the two tagging jets by identifying the jet pair with the highest invariant mass as the one produced by the scattering partons. This method correctly matches final state partons to the corresponding vector boson or scattering parton for at least 75% of the events under investigation. Events where the final state is produced through a triple vector boson emission are vetoed in order to select the EW VBS channel.
We simulate the irreducible QCD-induced l + l − + 4j sample and include its EFT dependence in the results. Another major background for this channel is Z+jets, which, however, is not included here due to the significant computational challenges posed by its Monte Carlo simulation. To prevent this choice from introducing a bias in the global analysis, the ZV+2j results will be analysed separately from the other channels, see Sec. 5.1. Up to this caveat, the semileptonic ZV+2j channel is found to be very competitive in terms of sensitivity to EFT effects, as it combines the presence of Z bosons in dominant diagrams with a large cross section. For this reason, it yields the strongest constraints on c HD , c (1) Hl and c HW B . It is also competitive with the inclusive WW on c Hq and c (1) ll , and with OSWW+2j on c HW . Distributions for interesting variables for some operators this process is sensitive to, after the selections, are shown in Figure 16.

Results
In this section we report the results of the likelihood scans with different configurations, corresponding to the statistics expected after an integrated luminosity of 100 fb −1 , collected by one given experiment at the LHC. Figure 2 shows −2∆ log L, profiled over the systematic nuisance parameter, as a function of individual Wilson coefficients, for each final state and for their combination (excluding ZV+2j). The 68% and 95% c.l. intervals obtained are reported in Figure 3 for the leptonic channels, and in Figure 4 for the semileptonic ZV+2j final state. In the latter two figures, the bands show the limits obtained with the baseline analysis, while the thin lines correspond to those obtained when neglecting all quadratic terms in Eq. (2.4). Note that the horizontal scale is not linear, in order to better visualise very different ranges in the image.

One-dimensional constraints
Optimal observables. As explained in 3.3, the likelihoods are constructed differently for each Wilson coefficient, by choosing optimal differential observables. For the 1D fits, the list of most sensitive observables is reported in Table 5.
Focusing on quadratic fits, for the four-quark operators, the most sensitive variables are found to be the transverse momentum of the first and second leading quarks p T,j 1 , p T,j 2 . This is expected, as these operators intuitively only affect the kinematics of the tagging jets. On the other hand, for the bosonic parameters c HW , c H , c W , that mainly enter the partonic scattering process, the variables that play the most important role are related to the kinematic properties of the leptonic final state. Transverse momentum variables related to the hadronic decay of the vector boson are the most sensitive to bosonic operators in the ZV+2j channel. Impact of individual processes. As expected, the inclusive WW gives the strongest constraints for most operators. Nevertheless, the VBS topology is found to be generally competitive. In the cases of four-quark and Higgs (c HW , c H ) operators, VBS is by construction the only constraining channel. In addition, as discussed in Sec. 4, OSWW+2j, ZZ+2j and WZ+2j are particularly competitive for c (1) Hl , c HD and c HW B . Globally, among the leptonic VBS processes, the most sensitive to EFT effects are SSWW+2j and OSWW+2j. This is mainly due to their cross-sections being larger than for WZ+2j and ZZ+2j, that are suppressed by the Z → branching ratio. Constraints obtained from the ZV+2j final state are also very competitive with the inclusive diboson ones, justifying the interest in this channel for a EFT analysis, and will deserve a more detailed study that includes also the backgrounds due to the production of a single vector boson plus jets.
Results by operator. In the combined analysis, the most constrained coefficients are the four-quark interactions c 1) qq , the parameter c W , that modifies TGC and QGC in a highly momentum-enhanced way, and the coefficients c ll . All these are bound to be below 0.15 at 68% c.l. for Λ = 1 TeV. This result is roughly consistent with the findings of Ref. [34] for the ZZ+2j case.
The four-quark operators differ among them by the SU(2) and flavour structures. In particular c   For each coefficient, the likelihood was built taking, for each channel, the distribution in the most constraining variable at 68% c.l. (see Table 5). Only the shown Wilson coefficient is varied at each time, and the others are set to 0. The sensitivity estimate for the OSWW+2j, WZ+2j, ZZ+2j and ZV+2j channels includes contributions from the respective QCD induced processes.
Linear+Quadratic 68% C.L. Linear+Quadratic 95% C.L.       give particularly momentum-enhanced signals in the longitudinally-polarized component and c W in the transverse one [70]. c (1) ll only enters via corrections to the EW input quantities. Its dominant effect is a rescaling of the overall cross sections, with nearly no shape modification in the distributions. The strongest bounds on this operator come again from diboson. c HD , c HW B also dominantly enter via input corrections and lead to qualitatively similar effects, but their impact in the inclusive WW and SSWW+2j is smaller compared to c (1) ll , resulting in weaker constraints. The coefficients c H and c HW only affect HVV couplings and, among the processes considered here, they are dominantly constrained by OSWW+2j, which is the one with the largest cross-section with the Higgs entering in the s-channel. These constraints are of course weaker than those imposed by available Higgs production and decay measurements.
Finally, the coefficient c Hl is significantly constrained in the ZZ+2j and WZ+2j VBS, which only give mild bounds.
Linear vs quadratic EFT parameterization. As we consider kinematic distributions extending to high energies, the validity of the EFT expansion cannot be guaranteed a priori over the entire parameter and phase spaces explored. Rather than implementing unitarisation or clipping procedures on the simulated events, that could introduce a dependence of the results on the specifics of these techniques, we provide a simple, qualitative assessment of the EFT validity by performing a comparison of the limits obtained retaining vs. neglecting quadratic SMEFT contributions (and omitting propagator corrections in both cases). Although exceptions are possible, a dominance of quadratic terms in the fit typically points to a breakdown of the EFT expansion, and indicates a potential sensitivity to neglected higher-dimensional operators. In Figs. 3, 4, one observes that quadratic terms significantly impact the combined results for less than a half of the operators, while, for the others, their inclusion has little consequence. Among the latter, c Hq , c Hl , c (1) ll are dominantly constrained through their linear contributions in all the VBS processes as well as in the inclusive WW. This is consistent, for instance, with the corresponding plots in Figs. 11, 13, 15, 16, which indicate that the linear contribution dominates in most bins already for c α /Λ 2 = 1 TeV −2 .
In the case of c Hl , c HD , c H , the sensitivity to quadratic terms varies between processes. However, it is very limited for the semileptonic channel and for the leptonic channels that dominate the combined constraints, i.e. WZ+2j for c (1) Hl (see Fig. 14) and OSWW+2j for c HD and c H (see Fig. 13). The case of c HW B is slightly different as several channels concur in constraining this operator. At quadratic level, the dominant constraint is WW. At linear level, this particular limit weakens significantly, but the effect is compensated by the constraints from ZZ+2j, WZ+2j and OSWW+2j, which become dominant and leave the final result nearly unchanged.
The four-quark Wilson coefficients c show the opposite behavior. All these operators are dominantly constrained by SSWW+2j, where they are measured in the p T,j i distributions. As shown in Fig. 12, for c α /Λ 2 = 1 TeV −2 , the interference contributions for all four-quark operators are negative and close in size to the quadratic terms. Intuitively, for Wilson coefficients sufficiently small the interference becomes largest and drives the fit. This threshold lies within the ballpark of the 95% c.l. sensitivity, and the fit happens to be such that the upper bounds on c are not. This effect is also responsible for the asymmetry of the constraints on these operators, seen in Figs. 2, 3: the bound for c α > 0 is weaker because it lies within a region where large cancellations take place between interference and quadratic terms. Indeed, neglecting the quadratic contribution spoils the cancellation and results in an improvement of this constraint, that is most visible for c (3,1) qq .
Finally, the coefficients that show the largest differences in the constraints derived with and without the inclusion of quadratic terms are c (1) qq , c (1,1) qq , c W , c HW and c (1) Hq . For all these, the difference is generally due to the quadratic contribution being significantly enhanced compared to the interference in the most relevant processes. This is can indeed be observed in the distribution plots shown in Appendix B (all the distributions can be found in ).
Impact of SMEFT corrections in propagators. As discussed in Sec. 2, SMEFT operators can contribute via both vertex and propagator corrections, due to the presence of EFT contributions to the total widths of the W, Z and Higgs bosons. Propagator corrections were omitted in the baseline results discussed so far, because currently available simulation tools only allow their consistent estimate at linear order in the Wilson coefficients. Therefore, we examine their impact within a purely linear setup: Figure 5 shows a comparison of the results obtained in individual fits retaining vs. neglecting propagator contributions. For consistency, the latter were added to both the EW and QCD-induced processes, and the two fits are performed using the same optimal variables in both cases. In fact we have verified that the inclusion of propagator effects does not alter the variables sensitivity ranking.
The width corrections are only sizeable in the phase space region where the intermediate boson is approximately on-shell. For this reason, contributions from δΓ H are very suppressed in all channels except OSWW+2j, and, in general, we observe that in each channel the most relevant propagator corrections are those associated to the reconstructed bosons.
The constraints that are most impacted by the inclusion of propagator corrections are those on c Hl . This is expected because these Wilson coefficients enter the total widths of both W and Z, and therefore enter all diagrams. In addition they give the largest contributions to both δΓ W and δΓ Z , see Eq. (3.2). 3 These corrections always partially cancel against the vertex contributions, leading to a large worsening of the constraints for all channels, that reaches a factor 2 in most cases and even a factor 7.5 for c (3) Hl in inclusive WW. It is important to underline that for such large cancellations to take place consistently across the entire fitted distribution, the vertex and propagator corrections need to have similar shapes. In the NWA, the latter amounts to an overall rescaling of the SM prediction in the resonant region. Therefore, cancellations can only happen for operators whose vertex corrections also "factor" into a SM rescaling when the intermediate boson is on-shell. This is indeed the case for c   the SM Higgs-mediated process: c HW induces a kinematic enhancement compared to the SM, and c HD , c HW B enter several non-Higgs diagrams.
In ZZ+2j and WZ+2j, the operators c (1) Hq , c (1) Hl , c HD and c HW B enter exclusively through δΓ Z . The constraints on c (1) Hq and c HW B change by over a factor of 2 in both processes, although in opposite directions. The constraints on c HD and c (1) Hl vary much less significantly, because they only give smaller contributions to δΓ Z , see Eq. (3.2).
Among the combined results, only those on c Hl , c (1) ll and C H are affected significantly by the introduction of propagator corrections, as most variations observed in individual channels either compensate each other, or are over-ruled by constraints from processes that are insensitive to propagator corrections (e.g. the constraint on c (1) Hq is still dominated by inclusive WW).
Impact of the QCD-induced sample. The QCD-induced processes are typically considered as a background to EW VBS analyses. Nevertheless, they are also generally modified in the presence of EFT operators. An interesting question is whether the dependence of the QCD background on the Wilson coefficients can have any significant impact on the constraints extracted from the statistical analysis. Here we compare the one-dimensional baseline results reported in Sec. 5.1 with those obtained fixing the QCD component to its SM shape. Figure 6 shows that for all operators and in all channels, the inclusion of the QCD EFT dependence never weakens the sensitivity: in some cases its impact is negligible, but in many others it leads to an improvement of the constraints by up to a factor of 2. In first approximation, this behavior can be understood considering that a worsening of the constraints could only occur if (partial) cancellations between EW and QCD SMEFT contributions took place in all bins of a distribution. However, for the vast majority of the operators and observables considered, this is not the case.
As a general rule, constraints for which the EW+QCD and EW results are very close, are numerically dominated by the EW-induced process alone. Among the constraints that show a significant improvement with the introduction of QCD, those on c (3) Hq and c (1) Hq in W ± W ∓ + 2j are strongly dominated by the pure QCD-induced contribution, while the others results from a non-trivial interplay between the EW and QCD bounds. These conclusions were checked against fits to the QCD-induced components only.

Individual vs. profiled bounds
When interpreting the EFT formalism as a low energy manifestation of a UV complete theory, it has been shown that its matching to the SMEFT basis typically introduces more than one non-zero Wilson coefficient [79,80]. The precise mapping of each subset of SMEFT operators to a UV complete model strictly depends on the model tested. In this context, upper bounds to the sensitivity can be estimated in the worst case scenario where all the operators under study are present simultaneously with non-zero Wilson coefficients.
We compare the results in Sec. 5.1 to those obtained allowing all 14 Wilson coefficients to float in the likelihood maximisation, and profiling over all of them except the one of interest. For each operator, the template analysis of the profiled constraint is carried out using, for each process, the same best variable obtained from individual constraints described in Section 5.1.
The comparison is shown in Figure 7. By definition, the profiled constraints are always equal to or worse than the individual ones. The differences observed between the two are quite heterogeneous, and they vary between a factor 1 (i.e. no difference) and 20. In particular, the constraints on c W , c qq , c (3,1) qq are nearly unaffected by the introduction of extra degrees of freedom, suggesting that these directions are well-resolved already within the individual processes that dominate the bounds, i.e. inclusive WW and SSWW+2j for the first and last two respectively. The largest deterioration in the constraints is observed for c (1) ll and c Hl . This can be easily traced back to a combination of the corresponding Wilson coefficients, approximately close to (c ll ), remaining nearly unconstrained in the fit. This is clearly visible in Figs. 8 and 9 and further discussed below.

Two-dimensional constraints
In this section we discuss constraints obtained allowing two operators to vary at the same time, fixing the remaining ones to zero. The analysis follows a strategy analogous to the one employed for the individual studies, as described in Sections 3.2 and 5.1. As above, EFT contributions to QCD-induced components of the VBS processes are accounted for in the fit,   Figure 7. Sensitivity of the combined analysis of VBS SSWW+2j, OSWW+2j, WZ+2j, ZZ+2j and diboson WW to the dimension-six operators considered, when the remaining Wilson coefficients are set to zero (green) or profiled away (grey). The QCD-induced EFT dependence was included where relevant.
whenever pertinent. The list of optimal observables employed for each channel and operator pair is provided in Tables 6-10 in Appendix A. Figures 8 and 9 show a subset of the likelihood scans obtained. The first figure illustrates the interplay between different processes for fixed operator pairs, while the second compares the combined sensitivity among all coefficient pairs containing c Hq and c Hl . The complete set of scans is available at the GitHub repository .
Interplay between measurements. For most operators, the profiled bounds are always dominated by a single process, that coincides with the leading constraint in the 1-dimensional fit. This is the case, for instance, of c HW , c H (OSWW+2j), c W , c Hq , c Hq , c Hl , c ll (WW) 4 and the 4-quark operators (SSWW+2j). Examples are shown in the (c W , c HW ) and (c HW , c HW B ) panels in Fig. 8, where the diboson and VBS constraints are nearly orthogonal to each other. The interplay between different channels is found to be most relevant for c HD and c HW B , see the corresponding panel in Fig. 8. All channels play a role in constraining these operators. The dominant bounds vary depending on the operator pair considered, although in most cases 4 The constraints on c Hl and c (1) ll are always dominated by diboson, except when the pair formed by these two operators is considered, as in this case an unconstrained direction emerges. See below.  Resolving degeneracies. Studying bi-dimensional likelihoods can provide insights about how potential degeneracies between operators are resolved in the fit. These mainly arise in two cases: among four-quark operators in VBS and between c (3) Hl and c (1) ll in all processes. All 4-quark operator pairs exhibit an unconstrained direction in the linear case, which is removed by the introduction of quadratic terms, see e.g. the (c Fig 8. Moreover, for all pairs, the combined constraint is dominated by the SSWW+2j channel alone, consistently with the discussion in Sec. 5.2. The degeneracy of singlet vs triplet SU (2) contractions (i.e. c ) is generally better resolved compared to the degeneracy between different flavor contractions, and leads to stronger projected constraints. This is again consistent with the 1-dimensional result (Sec. 5.1) and with the 2-dimensional fits shown in Fig. 9 (left panels), where the curves for c essentially unconstrained within the fit range. In the combined fit, this is not resolved by the inclusion of quadratic terms, but rather thanks to the unconstrained directions of the various processes having slightly different slopes. Specifically, the combined constraint is dominated by the interplay of inclusive WW and SSWW+2j, whose scattering amplitudes scale with (c ll ) respectively. This dependence follows from the interplay between corrections to vertices entering the W W production processes, that scale with ∆G F = 2c Linear vs quadratic EFT parameterization. Analogously to the individual limits above, we compare the combined results obtained with a linear or quadratic EFT parameterisation (neglecting all propagator contributions). In Fig. 8, these are indicated respectively by dashed and solid black lines. We find that the impact of quadratic terms on our fit is generally sizeable for all the processes considered.
Examining the bounds obtained profiling over one of the two parameters in each 2D fit, the behavior observed in the 2D fits is qualitatively consistent with the one observed in the individual case. The Wilson coefficients whose constraints vary the most depending on whether or not the quadratics are retained are c W and c (1) Hq . The 4-quark operators are most sensitive to the quadratic terms in 2D fits where they are paired with one another, due to the presence of the unconstrained direction in the linear case. When paired with other coefficients, only c Hq , c (1) ll are the least sensitive to the introduction of quadratics, as their projected bounds remain essentially unchanged in all 2D planes. For the remaining parameters (c HW , c HD , c HW B , c H ) the bounds obtained with and without quadratics generally differ, but the size of the variation depends on the operator pair considered. At the bi-dimensional level, the introduction of quadratic terms generally leads to a deformation of the likelihood contours, such that in many cases one can identify either parameter space regions (sufficiently away from the SM point) that are 68% c.l. allowed in the quadratic fit but not in the linear one, and regions for which the opposite is true. Two examples are shown in the (c (3) qq , c (1) qq ) and (c HD , c HW B ) panels in Fig. 8.

Projected constraints for LHC Run III and HL-LHC
The LHC Run III is expected to deliver roughly 200 fb −1 after 3 years of activity. Combined with the Run I and Run II statistics, the data set will then amount to more than 300 fb −1 . The Run III will be followed by a long shut-down of the accelerator in order to prepare the machine for its high luminosity phase (HL-LHC). The instantaneous luminosity will be increased up to about 5 − 7.5 × 10 34 cm −2 s −1 , allowing the LHC to deliver approximately 3000 fb −1 in 10 years of data taking. This section presents a sensitivity projection to these two scenarios (keeping the centre-of-mass energy of the proton collisions at 13 TeV in both cases).
The projections are obtained simply by scaling the expected number of events computed for the individual constraints in Sec. 5.1 by a factor that accounts for the increase in luminosity. No scaling of the constraint on the 2% luminosity uncertainty is applied. Figure 10 shows the individual exclusion ranges at 95% c.l. expected for the VBS-only and for the VBS and inclusive WW combinations at 100, 300, 3000 fb −1 , including the relevant QCD-induced contributions and quadratic terms in the EFT predictions.
The projection study highlights that the VBS combination should be able to constrain at parton-level all the operators to less than [−1, 1] at 95% c.l. by the end of the HL-LHC, while the inclusion of the inclusive WW channel further improves these bounds down to the [−0.5, 0.5] level, reaching a few percent for the most constrained operators.

Summary and conclusions
We have presented a comprehensive study of the sensitivity of VBS and diboson measurements to dimension-six SMEFT operators. The full set of plots for the differential distributions considered and likelihood scans performed are available at the dedicated GitHub repository github.com/MultibosonEFTStudies/D6EFTPaperPlots.
We estimated the expected limits on 14 dimension-six operators, that are representative of all classes of SMEFT corrections entering VBS and diboson. The baseline constraints are estimated for LHC measurement with √ s = 13 TeV and with an integrated luminosity of 100 fb −1 , working at parton-and tree-level, and including terms that are quadratic in the Wilson coefficients. We considered VBS processes as full 2 → 6 scatterings, including nonresonant diagrams in the event generation. This allowed us to perform a consistent template analysis of a significant number of differential distributions, for variables typically exploited in realistic LHC analyses. This study indicates that four-quark operators are mostly probed by the kinematics of the two tagging jets, while leptonic variables provide the best constraints on bosonic operators, see App. A. Interestingly, angular observables are often among the most Hq or c Hl . The VBS W + W − +2j, W ± Z +2j, ZZ + 2j channels include the respective QCD-induced processes. Only two Wilson coefficients are varied at a time, while the others are fixed to zero. Quadratic EFT contributions are included in all cases. In each panel, the operator on the y axis is fixed, while the one on the x axis is different for each curve (see legend).  Figure 10. Expected 95% c.l. constraints on individual SMEFT operators for an integrated luminosity of 100 fb −1 (grey), 300 fb −1 (red) and 3 ab −1 (blue). Constraints expected from the combination of the VBS channels SSWW+2j, OSWW+2j, WZ+2j, ZZ+2j are depicted as filled colored boxes while the combination of VBS and WW channels as hollow boxes with a colored border line. EFT effects in the QCD-induced processes are included whenever relevant. The result for VBS+WW is not indicated for operators that do not enter WW. In this case the limits from the VBS+WW combination would coincide with those from VBS only.
sensitive ones, indicating that the phenomenology of SMEFT effects in VBS processes is quite rich, and it goes beyond enhancements in transverse momenta or invariant mass distributions. Overall, we find that VBS processes are particularly sensitive to four-quark operators, that, individually, can be constrained even more strongly than operators affecting TGC and QGC. Operators correcting fermion-gauge interactions also give large effects to the VBS distributions, often with marked shape distortions. On the other hand, the sensitivity to operators modifying the Higgs-gauge couplings is poorer, due to Higgs contributions being suppressed by the selection cuts.
The constraints were derived for one and for two Wilson coefficient at a time, and also for each coefficient after profiling over the other 14. The comparison between individual and profiled constraints has shown degradation in sensitivity up to an order of magnitude. The results obtained in the 2D-scans highlight the complementarity between the different channels considered, as well as between VBS and diboson WW measurements. In most cases, combining all channels increases significantly the sensitivity to BSM effects. The impact of quadratic EFT contributions in the fit was also assessed, by comparing individual constraints obtained with and without retaining them in the EFT parameterization. We found that the inclusion of quadratic terms enhances considerably the constraining power of both the VBSonly and VBS + WW combinations for the dimension-six operators that induce the most energy-enhanced effects. On the other hand, for about a half of the operators considered, the linear contributions are dominant.
Limited to the linear analysis, we also investigated the impact of SMEFT corrections stemming from operator insertions in the propagators of intermediate W, Z and H states. We found that, although these contributions can have a strong impact for individual channels, in the combined analysis their introduction leads to a mild worsening of the constraints for three operators only, namely c Hl and c H . The inclusion of propagator corrections in quadratic fits is not possible at present, due to ambiguities in the distinction between conventional quadratic terms and double insertions.
Where relevant, our analysis kept into account EFT corrections to the QCD-induced background processes at order α 4 EW α 2 QCD . Interestingly, we find that the inclusion of EFT effects in the QCD-induced processes never weakens the constraints on individual Wilson coefficients. On the contrary, in several cases is enhances the sensitivity by up to a factor 2.
The sensitivity of the ZV+2j process with a semileptonic final state was studied separately, due to the lack of simulated Z+jets background events. The sensitivity of this process to BSM physics is largely unexplored in the literature, but was found to be competitive with that of inclusive diboson WW, at least at the level of individual constraints. Although this result should be taken as optimistic, due to the incomplete background treatment, it is certainly promising and worth further investigations in the future.
Finally, we derived projections for the individual constraints at the integrated luminosities expected for the end of the LHC Run III (300 fb −1 ) and of the HL-LHC (3 ab −1 ). In the latter scenario the VBS combination is expected to constrain all the operators in the (−1, 1) range at the 95% c.l. while the addition of the diboson WW channel improves the constraints to the (−0.5, 0.5) range.
Our study is a first step in the exploration of the potential of VBS measurements to constrain EFT effects. It can be improved and extended in several ways, for instance by including a larger set of dimension-six operators in the statistical analysis, or by fully accounting for the hadronisation and detector effects. A more refined analysis could also consider further (reducible) background sources, particularly for the semileptonic ZV+2j channel, where the Z+jets background was neglected altogether. One interesting question, that we leave for a future work, is whether the kinematic richness of VBS processes could help discriminating between different operators, for instance by exploiting angular information or by getting access to the polarisations of the weak bosons in final state.
Op.        Table 4, for various observables and Wilson coefficients in the SSWW+2j-EW process, at an integrated luminosity of 100 fb −1 . Solid lines show the total prediction for one Wilson coefficient at a time, with c α /Λ 2 = 0.01 (red), 0.1 (orange) or 1 TeV −2 (blue). The pure interference (quadratic) EFT component, normalized to c α /Λ 2 = 1 TeV −2 , is indicated with a purple (green) dashed line. For all distributions, the last bin comprises all the overflow events.  Figure 13.
Comparison of SM (filled histograms) and BSM (lines) expected number of events after the event selection in Table 4, for various observables and Wilson coefficients in the OSWW+2j-EW and OSWW+2j-QCD processes, at an integrated luminosity of 100 fb −1 . The SM distribution is represented as a stacked histogram, summing EW (light grey) and QCD (dark grey) components. Solid lines show the total prediction for one Wilson coefficient at a time, with c α /Λ 2 = 0.01 (red), 0.1 (orange) or 1 TeV −2 (blue). The pure interference (quadratic) EFT component, normalized to c α /Λ 2 = 1 TeV −2 , is indicated with a purple (green) dashed line. For all distributions, the last bin comprises all the overflow events.   Figure 14. Comparison of SM (filled histograms) and BSM (lines) expected number of events after the event selection in Table 4, for various observables and Wilson coefficients in the WZ+2j-EW and WZ+2j-QCD processes, at an integrated luminosity of 100 fb −1 . The SM distribution is represented as a stacked histogram, summing EW (light grey) and QCD (dark grey) components. Solid lines show the total prediction for one Wilson coefficient at a time, with c α /Λ 2 = 0.01 (red), 0.1 (orange) or   Figure 15. Comparison of SM (filled histograms) and BSM (lines) expected number of events after the event selection in Table 4, for various observables and Wilson coefficients in the ZZ+2j-EW and ZZ+2j-QCD processes, at an integrated luminosity of 100 fb −1 . The SM distribution is represented as a stacked histogram, summing EW (light grey) and QCD (dark grey) components. Solid lines show the total prediction for one Wilson coefficient at a time, with c α /Λ 2 = 0.01 (red), 0.   Table 4, for various observables and Wilson coefficients in the ZV+2j-EW and ZV+2j-QCD processes, at an integrated luminosity of 100 fb −1 . The SM distribution is represented as a stacked histogram, summing EW (light grey) and QCD (dark grey) components. Solid lines show the total prediction for one Wilson coefficient at a time, with c α /Λ 2 = 0.01 (red), 0.1 (orange) or 1 TeV −2 (blue). The pure interference (quadratic) EFT component, normalized to c α /Λ 2 = 1 TeV −2 , is indicated with a purple (green) dashed line. For all distributions, the last bin comprises all the overflow events.