Differential cross-section measurements for the electroweak production of dijets in association with a Z boson in proton–proton collisions at ATLAS

Differential cross-section measurements are presented for the electroweak production of two jets in association with a Z boson. These measurements are sensitive to the vector-boson fusion production mechanism and provide a fundamental test of the gauge structure of the Standard Model. The analysis is performed using proton–proton collision data collected by ATLAS at s=13TeV\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{s}=13\ \hbox {TeV}$$\end{document} and with an integrated luminosity of 139fb-1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$139\ \hbox {fb}^{-1}$$\end{document}. The differential cross-sections are measured in the Z→ℓ+ℓ-\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Z\rightarrow \ell ^+\ell ^-$$\end{document} decay channel (ℓ=e,μ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\ell =e,\mu $$\end{document}) as a function of four observables: the dijet invariant mass, the rapidity interval spanned by the two jets, the signed azimuthal angle between the two jets, and the transverse momentum of the dilepton pair. The data are corrected for the effects of detector inefficiency and resolution and are sufficiently precise to distinguish between different state-of-the-art theoretical predictions calculated using Powheg+Pythia8, Herwig7+Vbfnlo and Sherpa 2.2. The differential cross-sections are used to search for anomalous weak-boson self-interactions using a dimension-six effective field theory. The measurement of the signed azimuthal angle between the two jets is found to be particularly sensitive to the interference between the Standard Model and dimension-six scattering amplitudes and provides a direct test of charge-conjugation and parity invariance in the weak-boson self-interactions.


Introduction
Measurements that exploit the weak vector-boson scattering (VBS) and weak vector-boson fusion (VBF) processes have become increasingly prevalent at the Large Hadron Collider (LHC) in the last few years. In the Higgs sector, measurements of Higgs boson production via VBF have been used to determine the strength, charge-conjugation (C) and parity (P) properties of the Higgs boson's interactions with weak bosons [1-7]. These measurements have recently been augmented by the observation of the electroweak production of two jets in association with a weak-boson pair [8][9][10][11][12], which is extremely sensitive to the VBS production mechanism and provides a stringent test of the gauge structure of the Standard Model of particle physics (SM). In the search for physics beyond the SM, the VBF and VBS production mechanisms have been used to search for dark matter [13,14], heavyvector triplets [15], Higgs-boson pair production [16], and signatures of warped extra dimensions [17].
All of these measurements and searches rely on theoretical predictions to accurately model the electroweak processes that are sensitive to the VBF and VBS production mechanisms. Specifically, Monte Carlo (MC) event generators are used to optimise the event selection and to extract the electroweak signal from the dominant background, with the signal extraction typically performed using fits to kinematic spectra. However, it is known that the theoretical predictions from different event generators do not agree, both in the overall production rate [9] as well as in the kinematic properties of the final state [18]. Model-independent measurements that directly probe the kinematic properties of VBF and VBS are therefore crucial, to determine which event generators can be used reliably in physics analysis at the LHC experiments. This article presents differential cross-section measurements for the electroweak production of dijets in association with a Z boson (referred to as EW Z j j production). The EW Z j j process is defined by the t-channel exchange of a weak vector boson, as shown in Fig. 1a, b, and is very sensitive to the VBF production mechanism. Previous measurements of EW Z j j production by ATLAS [19,20] and CMS [21][22][23] have focused on measuring only an integrated fiducial crosssection in a VBF-enhanced topology. The analysis presented in this article measures differential cross-sections of EW Z j j production in the Z → + − decay channel ( = e, μ) and as a function of four observables; the transverse momentum of the dilepton pair ( p t, ), the dijet invariant mass (m j j ), the absolute rapidity 1 separation of the two jets (| y j j |), and the signed azimuthal angle between the two jets ( φ j j ). The φ j j variable is defined as φ j j = φ f − φ b , where the two highest transverse-momentum jets are ordered such that y f > y b [24]. Collectively, these four observables probe the important kinematic properties of the VBF and VBS production mechanisms. The measurements are performed using proton-proton collision data collected by the ATLAS experiment at a centre-of-mass energy of √ s = 13 TeV and with an integrated luminosity of 139 fb −1 .
The EW Z j j differential cross-section measurements presented here are sufficiently precise that they can be used to probe a diverse range of physical phenomena. First, under the assumption of no beyond-the-SM physics contributions to the EW Z j j process, the measurements can be used to distinguish between the SM EW Z j j predictions produced by different event generators or by different parameter choices within each event generators. In the short term, the measurements will therefore help determine which event generator predictions can be used reliably in analyses that seek 1 ATLAS uses a right-handed coordinate system with its origin at the nominal interaction point (IP) in the centre of the detector and the zaxis along the beam pipe. The x-axis points from the IP to the centre of the LHC ring, and the y-axis points upwards. Cylindrical coordinates (r, φ) are used in the transverse plane, φ being the azimuthal angle around the z-axis. The pseudorapidity is defined in terms of the polar angle θ as η ≡ − ln tan(θ/2), and is equal to the rapidity y ≡ 0.5 ln ((E + p z )/(E − p z )) in the relativistic limit. Angular distance is measured in units of R ≡ ( y) 2 + ( φ) 2 .
to exploit VBF and VBS at the LHC. In the longer term, the measurements will provide crucial input if the theoretical predictions are to be improved. Second, and more generally, the measurements provide a new avenue to search for signatures of physics beyond the SM. The differential cross-section as a function of φ j j , for example, is found to be particularly sensitive to anomalous weak-boson selfinteractions that arise from CP-even and CP-odd operators in a dimension-six effective field theory. This parity-odd observable has been proposed as a method to search for CPviolating effects in Higgs boson production [24], but has not yet been measured in a final state sensitive to anomalous weak-boson self-interactions.
The layout of the article is as follows. The ATLAS detector is briefly described in Sect. 2. The signal and background simulations used in the analysis are described in Sect. 3. The event reconstruction and selection are described in Sect. 4. The method used to extract the electroweak component is described in Sect. 5. This includes a data-driven constraint on the dominant background process in which the jets that are produced in association with the Z boson arise from the strong interaction (strong Z j j production) as shown in Fig. 1c, d. The corrections applied to remove the impact of detector resolution and inefficiency are described in Sect. 6. The experimental and theoretical systematic uncertainties are presented in Sect. 7. Finally, the differential cross-sections for EW Z j j production are presented in Sect. 8. Differential cross-sections for inclusive Z j j production are also presented in Sect. 8 for the signal and control regions used to extract the electroweak component. The EW Z j j differential cross-sections are used in Sect. 9 to search for anomalous weak-boson self-interactions. A brief summary of the analysis is given in Sect. 10.

ATLAS detector
The ATLAS detector [25] at the LHC covers nearly the entire solid angle around the collision point. It consists of an inner tracking detector surrounded by a thin superconducting solenoid, electromagnetic and hadronic calorimeters, and a muon spectrometer incorporating three large superconducting toroidal magnets.
The inner-detector system is immersed in a 2 T axial magnetic field and provides charged-particle tracking in the range |η| < 2.5. The high-granularity silicon pixel detector covers the vertex region and typically provides four measurements per track, the first hit normally being in the insertable B-layer (IBL) installed before the start of Run 2 [26,27]. The IBL is followed by the silicon microstrip tracker which usually provides eight measurements per track. These silicon detectors are complemented by the transition radiation tracker (TRT), which enables radially extended track recon- Fig. 1 Representative Feynman diagrams for EW Z j j production (a, b) and strong Z j j production (c, d). The electroweak Z j j process is defined by the t-channel exchange of a weak boson and at tree level is calculated at O(α 4 EW ) when including the decay of the Z boson. The strong Z j j process has no weak boson exchanged in the t-channel and at tree level is calculated at O(α 2 EW αs 2 ) when including the decay of the Z boson struction up to |η| = 2.0. The TRT also provides electron identification information based on the fraction of hits (typically 30 in total) above a higher energy-deposit threshold corresponding to transition radiation.
The calorimeter system covers the pseudorapidity range |η| < 4.9. Within the region |η| < 3.2, electromagnetic calorimetry is provided by barrel and endcap highgranularity lead/liquid-argon (LAr) calorimeters, with an additional thin LAr presampler covering |η| < 1.8, to correct for energy loss in material upstream of the calorimeters. Hadronic calorimetry is provided by the steel/scintillatortile calorimeter, segmented into three barrel structures within |η| < 1.7, and two copper/LAr hadronic endcap calorimeters. The solid angle coverage is completed with forward copper/LAr and tungsten/LAr calorimeter modules optimised for electromagnetic and hadronic measurements respectively.
The muon spectrometer comprises separate trigger and high-precision tracking chambers measuring the deflection of muons in a magnetic field generated by the superconducting air-core toroids. The field integral of the toroids ranges between 2.0 and 6.0 T m across most of the detector. A set of precision chambers covers the region |η| < 2.7 with three layers of monitored drift tubes, complemented by cathode-strip chambers in the forward region, where the background is highest. The muon trigger system covers the range |η| < 2.4 with resistive-plate chambers in the barrel, and thin-gap chambers in the endcap regions.
Interesting events are selected for further analysis by the level-one (L1) trigger system, which is implemented in custom hardware. The selections are further refined by algorithms implemented in software in the high-level trigger (HLT) [28]. The L1 trigger selects events from the 40 MHz bunch crossings at a rate below 100 kHz. The HLT further reduces the rate in order to write events to disk at about 1 kHz.

Dataset and Monte Carlo event simulation
The analysis is performed on proton-proton collision data at a centre-of-mass energy of √ s = 13 TeV. The data were recorded between 2015 and 2018 and correspond to an integrated luminosity of 139 fb −1 .
Monte Carlo event generators are used to simulate the signal and background events produced in the proton-proton collisions. These samples are used to optimise the analysis, evaluate systematic uncertainties, and correct the data for detector inefficiency and resolution. A summary of the event generators is presented in Table 1 and further details of each generator are given below.
Electroweak Z j j production was simulated using three MC event generators. The default EW Z j j sample was produced with Powheg-Box v1 [29-31] using the CT10nlo [32] parton distribution functions (PDF) and is accurate to nextto-leading order (NLO) in perturbative QCD. The sample was produced with the 'VBF approximation', which requires a t-channel colour-singlet exchange to remove overlap with diboson topologies [33]. The parton-level events were passed to Pythia 8.186 to add parton-showering, hadronisation and underlying-event activity, using the AZNLO [34] set of tuned parameters. The EvtGen program [35] was used for the properties of the bottom and charm hadron decays. This sample is referred to as Powheg+Py8 EW Z j j production.
The second EW Z j j sample was produced in the VBF approximation with Herwig7.1.5 [36,37]. The samples were produced at NLO accuracy in the strong coupling using Vbfnlo v3.0.0 [38] as the loop-amplitude provider. The MMHT2014LO PDF set [39] was used along with the default set of tuned parameters for parton showering, hadronisation and underlying event. EvtGen was used for the properties of the bottom and charm hadron decays. This sample is referred to as Herwig7+Vbfnlo EW Z j j production.
The third EW Z j j sample was produced in the VBF approximation with the Sherpa 2.2.1 event generator [40]. The samples were produced using leading-order (LO) matrix elements with up to two additional parton emissions. The NNPDF3.0nnlo PDFs [41] were used and the matrix elements were merged with the Sherpa parton shower using the MEPS@LO prescription [42]. Hadronisation and underlyingevent algorithms were used to construct the fully hadronic final state using the set of tuned parameters developed by the Sherpa authors. This sample is referred to as Sherpa EW Z j j production.
The dominant background arises from Z j j final states in which the two jets are produced from the strong interaction, as shown in Fig. 1c, d. This is referred to as the strong Z j j background and was simulated using three different MC event generators. Sherpa 2.2.1 was used to produce Z +nparton predictions (n = 0, 1, 2, 3, 4), at NLO accuracy for up to two partons in the final state and at LO accuracy for three or four partons in the final state, using the Comix [43] and OpenLoops [44,45] libraries. The different final-state topologies were merged into an inclusive sample using an improved CKKW matching procedure [42,46], which has been extended to NLO accuracy using the MEPS@NLO prescription [47]. The Sherpa prediction was produced using the NNPDF3.0nnlo PDFs and normalised to a next-to-nextto-leading-order (NNLO) prediction for inclusive Z -boson production [48]. The default set of tuned parameters in Sherpa was used for hadronisation and underlying-event activity. This sample is referred to as Sherpa strong Z j j production.
The second strong Z j j sample was produced using the MadGraph5_aMC@NLO generator [32] and is accurate to NLO in the strong coupling for up to two partons in the final state. The NNPDF2.3nlo PDF set [49] was used in the calculation. The MadGraph5_aMC@NLO generator was interfaced to Pythia 8.186 to provide parton showering, hadronisation and underlying-event activity, using the A14 set of tuned parameters. To remove overlap between the matrix element and the parton shower, the different jet multiplicities were merged using the FxFx prescription [50]. EvtGen was used for the properties of the bottom and charm hadron decays. The sample is normalised to the same NNLO prediction as for the Sherpa sample and is referred to as MG5_NLO+Py8 strong Z j j production.
The third strong Z j j sample was also produced with MadGraph5_aMC@NLO, but with the Z +n-parton matrixelements produced at LO accuracy for up to four partons in the final state. The NNPDF3.0lo PDFs were used in the calculation. The parton-level events were passed to Pythia 8.186 to provide parton-showering, hadronisation and underlyingevent activity, using the A14 set of tuned parameters [51]. To remove overlap between the matrix element and the parton shower, the CKKW-L merging procedure [52,53] was applied. EvtGen was used for properties of the bottom and charm hadron decays. The sample is normalised to the same NNLO prediction as for the Sherpa sample and is referred to as MG5+Py8 strong Z j j production.
Production of diboson (V V ) final states were simulated using Sherpa at NLO accuracy for up to one parton in the final state, and at LO accuracy for two or three partons in the final state. The NNPDF3.0nnlo PDF set was used in the calculation. The virtual corrections were taken from Open-Loops and the different topologies were merged using the MEPS@NLO algorithm. The default set of tuned parameters in Sherpa was used for hadronisation and underlying-event activity.
Backgrounds from events containing a single top quark or a top-antitop (tt) pair were estimated at NLO accuracy, using the hvq program [54] in Powheg-Box v2. The partonlevel events were passed to Pythia 8.230 to provide the parton showering, hadronisation and underlying-event activity using the A14 set of tuned parameters. EvtGen was used for the properties of the bottom and charm hadron decays. The NNPDF3.0nnlo PDF set was used and the h damp parameter in the Powheg-Box was set to 1.5 m top . The background from the W +jets final state was estimated using Sherpa, with the same set-up as for the Z +jets final state. The small contribution from triboson events (V V V production) was estimated using Sherpa at LO accuracy for up to one parton in the final state. The MEPS@LO prescription was used to merge the samples. The samples were produced using the NNPDF3.0nnlo PDF and the Sherpa authors' default parameterisation was used for hadronisation and underlyingevent activity.
The signal and background events were passed through the Geant4 [55] simulation of the ATLAS detector [56] and reconstructed using the same algorithms as used for the data (except for the Herwig7+Vbfnlo and MG5_NLO+Py8 samples, which were produced only at particle level). Differences in lepton trigger, reconstruction and isolation efficiencies between simulation and data are corrected on an event-by-event basis using p t -and ηdependent scale factors for each lepton [57,58]. The effect of multiple proton-proton interactions (pile-up) in the same or nearby bunch crossings is accounted for using inelastic proton-proton interactions generated by Pythia8 [59], with the A3 tune [60] and the NNPDF2.3LO PDF set [49]. These inelastic proton-proton interactions were added to the signal and background samples and weighted such that the distribution of the average number of proton-proton interactions in simulation matches that observed in the data.
An approximate detector-level prediction for MG5_NLO+Py8 is obtained by reweighting the strong Z j j simulation produced by MG5+Py8 such that the kinematic distributions match MG5_NLO+Py8 at particle level. This is referred to as MG5_NLO+Py8'. Similarly, an approximate detector-level prediction for Herwig7+Vbfnlo is obtained by reweighting the EW Z j j simulation produced by Powheg+Py8 to match Herwig7+Vbfnlo at particle level. This is referred to as Herwig7+Vbfnlo'.

Event reconstruction and selection
Events are required to pass unprescaled dilepton triggers with transverse momentum thresholds that depend on the lepton flavour and running periods. In 2015, the dielectron triggers retained events with two electron candidates that had p t > 12 GeV, whereas the dimuon triggers selected events with leading (subleading) muon candidates having p t > 18 (8) GeV. The transverse momentum thresholds for the lepton candidates were gradually increased during data taking, such that both electron candidates had p t > 24 GeV in 2018, whereas the leading muon threshold was increased to 22 GeV in the same running period.
Events are used in the analysis if they were recorded during stable beam conditions and if they satisfy detector and data-quality requirements [61]. The positions of the protonproton interactions are reconstructed using tracking information from the inner detector, with each associated vertex required to have at least two tracks with p t > 0.5 GeV. The primary hard-scatter vertex is defined as the one with the largest value of the sum of squared track transverse momenta.
Muons are identified by matching tracks reconstructed in the muon spectrometer to tracks reconstructed in the inner detector. Each muon is then required to satisfy the 'medium' identification criteria and the 'Gradient' isolation working point [57]. Muons are required to be associated with the primary hard-scatter vertex by satisfying |d 0 /σ d 0 | < 3 and |z 0 × sinθ | < 0.5 mm, where d 0 is the transverse impact parameter calculated with respect to the measured beam-line position, σ d 0 is its uncertainty, and z 0 is the longitudinal difference between the point at which d 0 is measured and the primary vertex. Reconstructed muons are used in the analysis if they have p t > 25 GeV and |η| < 2.4.
Electrons are reconstructed from topological clusters of energy deposited in the electromagnetic calorimeter that are matched to a reconstructed track [58]. They are calibrated using Z → ee data [62]. Each electron is required to satisfy the 'medium' likelihood identification criteria [58], as well as the same isolation working point as for muons. Electrons are required to be associated with the primary hard-scatter vertex by satisfying |d 0 /σ d 0 | < 5 and |z 0 × sinθ | < 0.5 mm. Reconstructed electrons are used in the analysis if they have p t > 25 GeV and |η| < 2.47, but excluding the transition region between the barrel and end-cap calorimeters (1.37 < |η| < 1.52).
Jets are reconstructed with the anti-k t algorithm [63,64] using a radius parameter of R = 0.4. The inputs to the algorithm are clusters of energy deposited in the electromagnetic and hadronic calorimeters. The jets are initially calibrated by applying energy-and pseudorapidity-dependent correction factors derived from simulation in the 'EM+JES' scheme [65], and then further calibrated using data-driven correction factors derived from the transverse momentum balance of jets in γ +jet, Z +jet and multijet topologies. Jets are used in the analysis if they have p t > 25 GeV and |y| < 4.4. As all highp t electrons pass the above requirements, jets are required to not overlap with a reconstructed Table 2 Observed and expected event yields in the dielectron and dimuon decay channels following the event selection described in Sect. 4. The first (second) uncertainty quoted for each generator is the experimental (theoretical) systematic uncertainty. The experimental systematic uncertainties are shown for each prediction. Theoretical uncertainties are calculated for all predictions except for MG5+Py8 strong Z j j, which is denoted 'N/A' in the Strong Z j j (Sherpa) 13 500 ± 600 ± 4500 15 100 ± 600 ± 5000 Strong Z j j (MG5+Py8) 13 140 ± 480 ± N/A 14 810 ± 540 ± N/A Strong Z j j (MG5_NLO+Py8') 8800 ± 300 ± 1000 10 000 ± 400 ± 1200 Other V V 45 ± 2 ± 2 4 5 ± 2 ± 2 tt, single top 92 ± 8 ± 6 9 8 ± 8 ± 6 W (→ ν)+jets, Z (→ τ τ )+jets Negligible Negligible electron (i.e. R( j, e) > 0.2). Jets with p t < 120 GeV and |η| < 2.4 are also required to be consistent with originating from the primary hard-scatter vertex using the 'medium' working point of the jet vertex tagger (JVT > 0.59) [66]. Following jet reconstruction, an additional quality requirement is placed on the events, by removing events containing jets that originate from noise bursts in the calorimeter. This removes 0.4% of the events in data.
Events are then selected if they have a topology consistent with EW Z j j production. A Z -boson candidate is reconstructed by requiring that each event contains exactly two charged leptons ( = e, μ) that are opposite in charge and of the same flavour. These leptons are required to be well separated from jets by imposing R( , j) > 0.4. The invariant mass and transverse momentum of the dilepton system is required to fulfil m ∈ (81, 101) GeV and p t, > 20 GeV. Events are required to contain two or more jets, with the leading and subleading jets satisfying p t > 85 GeV and p t > 80 GeV, respectively. The dijet system is then constructed from the two leading jets and is required to fulfil m j j > 1 TeV and | y j j | > 2.0. The Z boson is required to be centrally produced relative to the dijet system by imposing ξ Z < 1.0; the quantity ξ Z is defined as ξ Z = |y − 0.5(y j1 + y j2 )| / | y j j |, where y , y j1 and y j2 are the rapidities of the dilepton system, the leading jet, and the subleading jet, respectively. Finally, to reduce the impact of jets that originate from pile-up interactions and that survive the JVT selection criteria, the Z -boson candidate and the dijet system are required to be approximately balanced in transverse momentum, by requiring that p bal t < 0.15, where p bal t = | i p t,i | / i p t,i and the summation includes the dilepton system, the dijet system, and the highest transversemomentum additional jet reconstructed in the rapidity interval spanned by the dijet system.
The number of events in data that pass these selection requirements is shown in Table 2. The predicted event yield for each MC simulation is also presented. There is a large spread of EW Z j j event yields predicted by the different event generators. Furthermore, the predicted strong Z j j event yield also has significant uncertainties, with large theory uncertainties in each prediction and a large difference between the predictions of the different event generators. The contribution of the other processes amounts to about 3%.
The disagreement between data and simulation is not just observed in the total event yield. Figure 2 shows the data and predicted event yield as a function of m j j , | y j j |, p t, , and φ j j , with Sherpa used to model the strong Z j j process and Powheg+Py8 used for the EW Z j j process. The level of agreement between data and simulation depends on the kinematic properties of the event, with agreement at large m j j being particularly poor for this configuration of MC simulations.

Extraction of electroweak component
The poor agreement between data and simulation observed in Fig. 2 implies that the EW Z j j event yield cannot be extracted by simply subtracting the background simulations from the data. Furthermore, the level of mismodelling in the simulation changes when different strong Z j j simulations are used, as shown in Fig. 3 for the m j j and p t, distributions. A datadriven method is therefore used to constrain both the shape and normalisation of the strong Z j j background during the extraction of the EW Z j j event yield.
The data are split into four regions by imposing criteria on ξ Z as well as on the multiplicity of jets in the rapidity interval between the leading and subleading jets, N gap jets . These The m j j spectrum is shown starting from 250 GeV, and hence includes more events than the other plots that use the default m j j > 1000 GeV criterion Fig. 3 Ratio of Monte Carlo prediction to data for different physics processes and generators for the m j j and p t, distributions, following the event selection described in Sect. 4. The data contain all processes that pass the event selection and the ratio demonstrates the contribution to the observed event yield that is predicted by each MC generator. The m j j distribution extends down to 250 GeV and hence includes a larger phase space than the p t, distribution, which requires m j j > 1000 GeV. Only statistical uncertainties are shown. The prediction labelled MG5_NLO+Py8' for the strong Z j j prediction is obtained by a particle-level reweighting of the strong Z j j simulation provided by MG5+Py8. The EW Z j j prediction labelled Herwig7+Vbfnlo' is also obtained by a particle-level reweighting of the EW Z j j simulation provided by Powheg+Py8 Fig. 4 Definition of the signal region (SR) and control regions (CRa, CRb, CRc) used in the extraction of the electroweak component two variables are chosen because they are almost uncorrelated for both the strong and EW Z j j processes, with calculated correlation coefficients ranging from −0.04 to +0.02 depending on the event generator and process. Approximately 80% of the EW Z j j events are predicted to fall into the EW-enhanced signal region (SR) defined by N gap jets = 0 and ξ Z < 0.5. The remaining three regions define EWsuppressed control regions (CR), which can be used to constrain the dominant background from strong Z j j production. These regions are labelled as CRa (N gap jets ≥ 1, ξ Z < 0.5), CRb (N gap jets ≥ 1, ξ Z > 0.5) and CRc (N gap jets = 0, ξ Z > 0.5) and are depicted in Fig. 4. All analysis decisions and optimisations were performed with the signal region blinded, to avoid any unintended biases.
The EW Z j j event yield is measured in the EW-enhanced SR using a binned maximum-likelihood fit [67,68]. The log likelihood is defined according to where r is an index corresponding to the region r ∈ {CRa, CRb, CRc, SR}, i is the bin of the kinematic observable, N data ri is the observed event yield and ν ri (θ ) is the prediction that is dependent on the s sources of experimental systematic uncertainty that are each constrained by nuisance parameters θ = (θ 1 , . . . , θ s ). 2 The fitted number of events in each region and in each bin of a distribution is given by (1) 2 The dependence of the prediction on the systematic uncertainties is given by ν ri (θ) = ν MC ri s (1 + λ ris θ s ), where s is an index for the uncertainty source, θ s is the associated nuisance parameter and λ ris is the fractional uncertainty amplitude for bin i in region r .
where μ i is the EW Z j j signal strength of bin i, ν EW,MC ri and ν other,MC ri are the MC predictions of EW Z j j and contributions from other processes (diboson, tt and single top), respectively. The strong Z j j prediction is constrained using the different EW-suppressed control regions according to (2) Here, the b L,i and b H,i are sets of bin-dependent factors that apply to the ξ Z < 0.5 and ξ Z > 0.5 regions, respectively. These factors are primarily constrained in CRa and CRb, where they adjust the predicted simulated strong Z j j event yields and bring the total predicted yield (v ri of Eq. 1) into better agreement with data. The f (x i ) is a two-parameter function of the observable that is being measured and is evaluated at the centre of each bin. This function provides a residual correction to the constrained strong Z j j yield to account for the extrapolation from CRa (N gap jets ≥ 1) to the SR (N gap jets = 0) and is primarily constrained by CRb and CRc. The function is taken to be a first-order polynomial.
The free parameters in the binned maximum-likelihood fit are therefore the signal strengths μ i , the two parameters of the function f (x i ), and the b L,i and b H,i corrections to the strong Z j j process. In total, this amounts to 3 N bins + 2 parameters that are constrained using 4 N bins measurements in data, where N bins is the number of bins measured for a specific observable (m j j , | y j j |, p t, and φ j j ).
The pre-fit and post-fit agreement between data and simulation is shown in Fig. 5 as a function of m j j in the signal and control regions. Two separate fits are shown, one using the Sherpa strong Z j j prediction (top row) and one using the MG5_NLO+Py8' prediction (bottom row). These simulations initially have very different mismodelling as a function of m j j , but produce very good agreement with the data following the fitting procedure. The overall scaling factor applied to the strong Z j j prediction from MG5_NLO+Py8' in the signal region is 0.93 at low m j j rising to 2.2 at high m j j . For Sherpa, the corresponding scaling factors are 0.86 at low m j j and 0.26 at high m j j . The pre-fit systematic uncertainties shown on the plots are derived as outlined in Sect. 7.
Since there is no a priori reason to prefer any strong Z j j generator over another, the EW Z j j component is extracted three times, once using the Sherpa strong Z j j prediction, once using the MG5_NLO+Py8' strong Z j j prediction, and once using the MG5+Py8 strong Z j j prediction. The final electroweak signal yield in each bin of the differential distribution is taken to be the midpoint of the envelope of yields obtained using the three different strong Z j j event genera- 5 Comparison between data and prediction before (left) and after (right) the fit using strong Z j j estimates based on Sherpa (top) and MG5_NLO+Py8' (bottom) in bins of m j j in the different control and signal regions. The MG5_NLO+Py8' prediction is obtained by a particlelevel reweighting of the strong Z j j simulation provided by MG5+Py8. The m j j bin edges are defined by (1.0, 1.5, 2.25, 3.0, 4.5, 7.5) TeV tors. The envelope itself is used to define a systematic uncertainty as outlined in Sect. 7.
The constraints on the strong Z j j simulation in Eq. 2 are evaluated independently for each of the measured differential distributions (m j j , | y j j |, p t, and φ j j ). This results in slightly different total EW Z j j and strong Z j j event yields when summed across each differential spectrum. To ensure consistency between the distributions, an additional constraint is applied in the likelihood to ensure that the same integrated strong Z j j yield is obtained for each distribution, i.e.
whereν strong SR,m j j is the event yield obtained by integrating the constrained strong Z j j template for the m j j distribution in the SR .
The electroweak extraction methodology is validated in four ways. First, a variation of the likelihood method is implemented by switching the control regions used to define the strong Z j j simulation as defined in Eq. 2, such that the b i factors are constrained in CRs at high ξ Z and the f (x i ) func-tion is then defined to correct for non-closure when transferring these corrections to low ξ Z . Second, the constraint on the strong Z j j background includes a function ( f (x i )) that is taken to be a first-order polynomial by default. This choice is validated by changing the function to a secondorder polynomial. Third, the constraint applied to the integrated strong Z j j event yield (Eq. 3) is removed. Finally, a simpler 'sequential' method is used to extract the EW Z j j event yields. In this approach, the data-driven correction to the strong Z j j is derived in CRa (assuming the SM prediction for the electroweak process in this region) and directly applied to the strong Z j j simulation in the SR. A transfer factor to account for mismodelling between the SR and CRa is evaluated at low m j j (250 ≤ m j j < 500 GeV). Nonclosure of the sequential method is evaluated in CRc using corrections to the strong Z j j process derived in CRb; this non-closure is used as a systematic uncertainty in the sequential method. The extracted electroweak event yields obtained with these four variations are found to be in good agreement with the nominal results and are presented in Appendix A.

Correction for detector effects
Particle-level differential cross-sections are produced by correcting the inclusive Z j j and EW Z j j event yields in each bin for the effects of detector inefficiency and resolution. The EW Z j j event yields are extracted in the signal region using the method outlined in the previous section. The inclusive Z j j event yields are obtained by subtracting, from the data, the small number of events predicted by simulation for processes that do not contain a Z boson and two jets in the final state (tt, single-top, V V → Z j j, and W +jets production). For both inclusive and EW Z j j production, the event yields in the e + e − and μ + μ − decay channels are added together and unfolded in a single step. The particle level is defined using final-state stable particles with mean lifetime satisfying cτ > 10 mm. To reduce model-dependent extrapolations across kinematic phase space, the particle-level event selection is defined to be as close as possible to the detector-level event selection defined in Sect. 4. Leptons are defined at the 'dressed' level, as the four-momentum combination of a prompt electron or muon (that do not originate from the decay of a hadron) and all nearby prompt photons within R < 0.1. Leptons are required to have p t > 25 GeV and have the same acceptance requirement as used at the analysis level, i.e. muons satisfy |η| < 2.4 and electrons satisfy |η| < 2.47 (but exclude the region 1.37 < |η| < 1.52). Jets are reconstructed using the anti-k t algorithm using all final-state stable particles as input, except those that are part of a dressed-lepton object. Jets are required to have p t > 25 GeV and |y| < 4.4. Using these jets and leptons, events are then selected in a VBF topology using requirements identical to those imposed at detector level. The EW Z j j differential cross-sections are measured in the SR, whereas inclusive Z j j differential cross-sections are measured in the SR and the three CRs. The VBF topology, SR and the three CRs are defined in Table 3.
Each distribution is unfolded separately using the iterative Bayesian method proposed by D'Agostini [69,70] with two iterations. This procedure uses MC simulations to (i) correct for events that pass the detector-level selection but not the particle-level selection, (ii) invert the migration between bins of the differential distribution, and (iii) correct for events that pass the particle-level selection but not the detectorlevel selection. For the EW Z j j differential cross-section measurements, the Powheg+Py8 EW Z j j simulation is used to define the corrections and the response matrices. For the inclusive Z j j differential cross-section measurements, all sources of Z j j production are part of the measurement and the unfolding is carried out using the cross-section weighted sum of the Powheg+Py8 EW Z j j simulation, the Sherpa strong Z j j simulation, and the Sherpa diboson samples that contain a leptonically decaying Z boson produced in association with a hadronically decaying weak boson.
Statistical uncertainties in the data are propagated through the unfolding procedure using the bootstrap method [71] with 1000 pseudo-experiments. For the EW Z j j measurements, the electroweak extraction is repeated for each pseudoexperiment after fluctuating the event yields, in each bin of the signal and control regions, using a Poisson distribution. For the inclusive Z j j measurements, the backgroundsubtracted event yields are fluctuated using a Gaussian distribution centred on the data-minus-background value and with a width given by the data statistical uncertainty. The statistical uncertainties in the MC simulation are propagated through the unfolding procedure in a similar fashion, by fluctuating each bin of the response matrix using a Gaussian distribution. The unfolding is repeated with the modified distributions (or response matrices) created for each pseudo-experiment. The final statistical uncertainties in the measurement are taken to be the standard deviation of the unfolded values obtained from the ensemble of pseudo-experiments.

Systematic uncertainties
Experimental systematic uncertainties Experimental systematic uncertainties arise from jet reconstruction, lepton reconstruction, the pile-up of multiple proton-proton interactions, and the luminosity determination. These uncertainties affect the normalisation and shape of the background simulations used in the extraction of the EW Z j j process, as well as the MC simulations used to unfold the EW Z j j and inclusive Z j j event yields. For the extraction of the electroweak signal, each source of experimental uncertainty is included as a Gaussian-constrained nuisance parameter in the likelihood, as outlined in Sect. 5. For the unfolding, each source of uncertainty is propagated to the MC simulations and the change in the unfolded event yield is taken as the systematic uncertainty.
The luminosity is measured to an accuracy of 1.7% using van der Meer beam separation scans, as outlined in Refs. [72,73]. Uncertainties in the modelling of pile-up interactions are estimated by repeating the analysis after varying the average number of pile-up interactions in the simulation. This variation accounts for the uncertainty in the ratio of the predicted and measured inelastic cross-sections within the ATLAS fiducial volume [74].
A variation in the pile-up reweighting of simulated events (referred to as pile-up uncertainty) is included to account for the uncertainty in the ratio of the predicted and measured inelastic cross-sections.
The lepton trigger, reconstruction and isolation efficiencies in simulation are corrected using scale factors derived from data, as outlined in Sect. 3. Systematic uncertainties associated with this procedure are estimated by varying these scale factors according to their associated uncertainties [57,58]. In addition, uncertainties due to differences between data and simulation in the reconstructed lepton momentum [57,62] are estimated by scaling and smearing the lepton momentum in the simulation. The overall impact on the differential cross-section measurement from systematic uncertainties associated with leptons is typically 1%, but rises to 2% at the highest dilepton transverse momentum.
The uncertainties associated with jet energy scale and jet energy resolution have a larger impact on the analysis. As discussed in Sect. 4, the jets are calibrated in data using a combination of MC-based and data-driven correction factors. The uncertainty in the measurement due to these corrections is estimated by scaling and smearing the jet fourmomentum in the simulation by one standard deviation in the associated uncertainties of the calibration procedure [65]. The impact on the differential cross-section measurements is between 5% at low m j j or p t, , but more than 10% for m j j > 4 TeV. An additional uncertainty arises from the use of the jet vertex tagger, which suppresses jets arising from pile-up interactions but is not fully efficient for jets produced in the hard scatter. Uncertainties arising from imperfect modelling of the JVT efficiency are estimated by varying the JVT requirement [66] and result in an uncertainty of about 1%, which is anti-correlated between the N gap jets = 0 and N gap jets ≥ 1 regions.
Theoretical uncertainties in the electroweak signal extraction Theoretical uncertainties associated with the modelling of the signal and background processes can impact the extraction of the electroweak signal yield. The impact of each source of theory uncertainty on the extracted signal yield is evaluated by repeating the electroweak extraction procedure (outlined in Sect. 5) after varying the input MC event generator templates in the SR and the CRs. The variation in the extracted signal yield is then propagated through the unfolding procedure.
Theoretical uncertainties associated with the modelling of the strong Z j j process are the dominant uncertainties in the extraction of the electroweak signal yield. Three sources of uncertainty in the strong Z j j modelling are investigated, arising from (i) the choice of event generator, (ii) the renormalisation and factorisation scale dependence in the strong Z j j calculations, and (iii) the parton distribution functions. The systematic uncertainty associated with the choice of event generator is defined by the envelope of electroweak event yields extracted using the Sherpa, MG5_NLO+Py8' and MG5+Py8 strong Z j j simulations (the default electroweak event yield defined as the midpoint of this envelope, as discussed in Sect. 5). The uncertainty associated with the choice of renormalisation and factorisation scales is assessed by repeating the analysis using new strong Z j j templates for Sherpa in which the renormalisation (μ R ) and factorisation (μ F ) scales have been varied independently by factors of 0.5 and 2.0. Six variations are considered for each generator corresponding to (μ R , μ F ) = (0.5, 1.0), (2.0,1.0), (1.0, 0.5), (1.0, 2.0), (0.5,0.5) and (2.0,2.0). For each variation, the change in the extracted EW event yield relative to that obtained with the default Sherpa strong Z j j sample is evaluated, and the envelope of the variations is then taken to be the relative uncertainty in the extracted electroweak yields. Finally, the impact of uncertainties associated with the parton distribution functions is estimated using the Sherpa generator, by reweighting the nominal strong Z j j sample to reproduce the variations of the NNPDF3.0nnlo PDF set (including the associated αsvariations) and repeating the full analysis chain for each variation. The systematic uncertainty in the extracted EW signal yields due to PDFs is then taken as the RMS of signal yields extracted from the PDF set variations. Of the three sources of uncertainty associated with modelling strong Z j j production, the choice of event generator has the largest impact on the extracted electroweak yields.
Theoretical uncertainties associated with the modelling of the EW Z j j process have a much smaller impact on the extraction of the electroweak component because, for each bin of a measured distribution, the only theoretical input is the relative event yields in the SR and CRs. The theoretical uncertainty due to the mismodelling of the EW Z j j process is determined by repeating the analysis after reweighting the default Powheg-Box EW Z j j simulation such that it matches the prediction of the Herwig7+Vbfnlo EW Z j j simulation at particle level. The change in extracted EW event yield with respect to the nominal event yield extracted with Powheg+Py8 is taken as a symmetric uncertainty. The signal-modelling dependence is further validated using the leading-order Sherpa EW Z j j simulation to extract the electroweak event yield and the results are found to be consistent and within the assigned uncertainty due to electroweak Z j j modelling. Systematic uncertainties associated with the parton distribution functions used in the matrix-element calculation are investigated, by applying the NNPDF3.0nnlo PDF set variations to the Sherpa EW Z j j simulation, and found to have a much smaller impact than the choice of event generator. Variations of renomalisation and factorisation scales in the matrix-element calculations are also found to have a negligible impact on the final result. The total systematic uncertainty associated with the signal modelling is typically between 2-3%.
The electroweak extraction methodology assumes that there is no interference between the EW Z j j process and the strong Z j j process. The size of the interference contribution relative to the electroweak signal process is estimated at particle level using MadGraph5 as a function of the measured kinematic variables in the SR and CRs. The uncertainty associated with the interference is then defined as the change in the extracted electroweak yield induced by reweighting the default Powheg+Py8 EW Z j j sample such that it contains the interference contribution, and is taken to be symmetric. This source of uncertainty is typically a factor of five smaller than the uncertainty associated with the modelling of the strong Z j j process.
Uncertainties in the unfolding procedure Uncertainties associated with the unfolding procedure are estimated in two ways. First, the data are unfolded using a different simulation and the deviation from the nominal result is taken as a systematic uncertainty. For the EW Z j j differential cross-section measurements, the Sherpa EW Z j j simulation is used in place of the Powheg+Py8 EW Z j j simulation. For the inclusive Z j j differential cross-section measurements, the MG5+Py8 strong Z j j simulation is used in place of the Sherpa strong Z j j simulation. Second, a data-driven closure test is performed separately for each observable, to assess the potential bias in the unfolding method. In this approach, the particle-level distribution is reweighted such that it provides a better description of the data at detector level. The reweighted detector-level prediction is then unfolded using the response matrix and other corrections derived from nominal (unweighted) Powheg+Py8 EW Z j j simulation. The systematic uncertainty associated with the unfolding method is defined as the difference between the unfolded spectrum and the reweighted particle-level prediction; it is taken to be a symmetric uncertainty.

Summary of systematic uncertainties
The final uncertainties in the differential cross-section measurements of EW Z j j production and inclusive Z j j production are shown in Fig. 6. For the inclusive Z j j measurements, the jet energy scale and jet energy resolution uncertainties dominate. However, for the EW Z j j measurements the uncertainties associated with the modelling of the strong Z j j process dominate.

Results
The differential cross-sections for EW Z j j production as a function of m j j , | y j j |, p t, , and φ j j are shown in Fig. 7 and are compared with theoretical predictions produced by Herwig7+Vbfnlo , Powheg+Py8 and Sherpa. The set-up of the theoretical predictions is discussed in Sect. 3. The effects of scale uncertainties on the Her-wig7+Vbfnlo prediction are estimated by independently varying the scale used in the matrix-element calculation and the scale associated with the parton shower by factors of 0.5 or 2.0. The effects of scale uncertainties on the Sherpa prediction are estimated by varying the renormalisation and factorisation scales used in the matrix-element calculation independently by a factor of 0.5 or 2.0. The effects of scale uncertainties on the Powheg+Py8 prediction are evaluated by independently varying the renormalisation, factorisation and resummation scales by factors of 0.5 or 2.0. Additional uncertainties on the Powheg+Py8 prediction associated with the parton-shower and underlying-event parameters in Pythia8 are evaluated using the AZNLO eigentune variations [34]. PDF uncertainties on the EW Z j j predictions are estimated by reweighting the nominal sample to reproduce the 100 variations of the NNPDF3.0nnlo PDF sets and taking the RMS of these variations; the impact of PDFrelated uncertainties on the EW Z j j predictions are found to be much smaller than the impact of scale uncertainties.
In general, the Herwig7+Vbfnlo prediction is found to be in reasonable agreement with the data for all measured distributions. The Powheg+Py8 prediction is found Fig. 6 Fractional uncertainty in the inclusive Z j j measurement (top) and the EW Z j j measurement (bottom) as a function of m j j (left) and p t, (right). Uncertainty sources are grouped in categories that are added in quadrature (denoted ⊕) to give the total uncertainty. The 'EW Z j j model' component includes the uncertainty on the EW Z j j prediction and the impact of interference between the strong Z j j and EW Z j j processes. The 'strong Z j j model' uncertainty is dominated by the choice of generator used for the strong Z j j prediction, but also includes the impact of renormalisation/factorisation scale variations and PDF set variations to overestimate the EW Z j j cross-section at high m j j , high | y j j |, and intermediate p t, . Furthermore, the central value of the Powheg+Py8 prediction often does not agree with the Herwig7+Vbfnlo prediction, within the assigned theoretical uncertainties. A similar discrepancy between theoretical predictions was noted for EW V V j j processes in Ref. [18] and was attributed to the set-up of the parton shower when matched to the matrix-element calculations. The Sherpa prediction significantly underestimates the measured differential cross-sections, due to a non-optimal setting of the colour flow [18]. However, despite the offset in normalisation, the shape of the measured distributions is reasonably well produced by Sherpa. Under the assumption that there are no new physics contributions to the EW Z j j process, the measurements presented in this article therefore constrain the choice of theoretical predictions that should be used for signal modelling in future measurements that exploit weak-boson fusion or weak-boson scattering. In particular, the EW Z j j differential cross-section measurements can be used to determine the optimal parameter choices for each event generator, and poor parameter choices can be ruled out entirely.
A fiducial cross-section for EW Z j j production is calculated, by integrating the differential cross-section as a function of m j j , and found to be σ EW = 37.4 ± 3.5 (stat) ± 5.5 (syst) fb. This is in excellent agreement with the theoretical prediction from Herwig7+Vbfnlo , which is 39.5 ± 3.4 (scale) ± 1.2 (PDF) fb. The unfolded data are shown as black points, with the statistical uncertainty represented by an error bar and the total uncertainty represented as a grey band. The data are compared with theoretical predictions produced by Herwig7+Vbfnlo (red points), Powheg+Py8 (blue points) and Sherpa 2.2.1 (orange points). Uncertainty bands are shown for the three theoretical predictions. Each theory prediction is slightly offset from the bin center to avoid overlap Differential cross-sections for inclusive Z j j production as a function of m j j , | y j j |, p t, and φ j j are also measured in the signal and control regions that are used to extract the electroweak component. These measurements can be used to re-evaluate the electroweak contribution in the future, when new theoretical predictions for the strong Z j j background presumably will become available. The differential cross-sections for inclusive Z j j production measured in the SR as a function of m j j , | y j j |, p t, , and φ j j are shown in Fig. 8. The differential cross-sections measured in CRa for inclusive Z j j production as a function of m j j and p t, are shown in Fig. 9. The data are compared with the strong Z j j predictions provided by Sherpaand MG5_NLO+Py8, augmented with the EW Z j j contribution predicted by Herwig7+Vbfnlo and the V Z contribution predicted by Sherpa. The effects of scale uncertainties on the strong Z j j predictions dominate the overall uncertainty in each prediction and are estimated by independently varying the renormalisation and factorisation scales by factors of 0.5 and 2.0 (with six variations considered for each generator). PDF uncertainties on the strong Z j j predictions are estimated using the variations of the NNPDF PDF sets. The total uncertainty on the strong Z j j predictions is taken to be the envelope of the scale variations added in quadrature with the PDF uncertainty. Overall, the data is best described when using the MG5_NLO+Py8 prediction for strong Z j j production.
The unfolded differential cross-sections for EW Z j j production and inclusive Z j j production are documented in tabular form in Appendix B. The data are also provided in the HEPDATA repository [75] and a Rivet analysis routine is provided [76,77]. The unfolded data are shown as black points, with the statistical uncertainty represented by an error bar and the total uncertainty represented as a grey band. The data are compared with theoretical predictions constructed from different strong Z j j predictions provided by Sherpa (green) and MG5_NLO+Py8 (blue). Uncertainty bands are shown for the two theoretical predictions. Each theory prediction is slightly offset from the bin center to avoid overlap

Constraints on anomalous weak-boson self-interactions
In this section, the measured EW Z j j differential crosssections are used to constrain extensions to the SM that produce anomalous weak-boson self-interactions. The anomalous interactions are introduced using an effective field theory (EFT), for which the effective Lagrangian is given by where L SM is the SM Lagrangian, the O i are dimension-six operators in the Warsaw basis [78], and the c i / 2 are Wilson coefficients that describe the strength of the anomalous interactions induced by those operators. Constraints are placed on two CP-even operators (O W , O H W B ) and two CP-odd operators (Õ W ,Õ H W B ), which are known to produce anomalous W W Z interactions. Theoretical predictions are constructed for the EW Z j j process using the effective Lagrangian in Eq. 4. The amplitude for the EW Z j j process is split into a SM part, M SM , and a dimension-six part, M d6 , which contains the anomalous interactions. The differential cross-section or squared amplitude then has three contributions Fig. 9 Differential cross-sections measured in CRa for inclusive Z j j production as a function of m j j (left) and p t, (right), where CRa is defined by N gap jets ≥ 1 and ξ Z < 0.5. The unfolded data are shown as black points, with the statistical uncertainty represented by an error bar and the total uncertainty represented as a grey band. The data are compared with theoretical predictions constructed from different strong Z j j predictions provided by Sherpa (green) and MG5_NLO+Py8 (blue). Uncertainty bands are shown for the two theoretical predictions. Each theory prediction is slightly offset from the bin center to avoid overlap namely a pure SM term |M SM | 2 , a pure dimension-six term |M d6 | 2 , and a term that contains the interference between the SM and dimension-six amplitudes, 2 Re(M * SM M d6 ). The constraints on the dimension-six operators presented in this section are derived both with and without the pure dimensionsix terms included in the theoretical prediction. This tests whether the results are robust against missing dimensioneight operators in the EFT expansion.
The pure-SM contribution to the EW Z j j differential cross-sections in Eq. 5 is taken to be the prediction from Her-wig7+Vbfnlo . The contributions arising from the interference and pure dimension-six terms are generated at leading order in perturbative QCD using MadGraph5+Pythia8, with the interactions from the dimension-six operators provided by the SMEFTSim package [79]. The A14 set of tuned parameters is used for parton showering, hadronisation and multiple parton scattering. To account for missing higherorder QCD corrections, the interference and pure dimensionsix contributions are scaled using a bin-dependent K -factor, which is defined by the ratio of pure-SM EW Z j j differential cross-sections predicted by Herwig7+Vbfnlo and Mad-Graph5+Pythia8 in each bin.
The impact of the interference and pure dimension-six contributions to the EW Z j j differential cross-sections is shown relative to the pure SM contribution in Fig. 10. The Wilson coefficients were chosen to be c W / 2 = 0.2 TeV −2 , c W / 2 = 0.2 TeV −2 , c H W B / 2 = 1.8 TeV −2 and c H W B / 2 = 1.8 TeV −2 . For the CP-even O W operator, the highp t, region is particularly sensitive to the anoma-lous interactions, a feature that was seen in previous studies for EW V j j production [23,80]. The pure dimension-six contributions to the cross-section dominate in this region. The φ j j observable is also found to be very sensitive to the anomalous interactions induced by the O W operator, but in this observable the interference contribution dominates. For the CP-even O H W B operator, the interference contribution dominates in all distributions, with the φ j j observable showing the largest kinematic dependences. For the CP-odd operators, the interference contribution is zero in the parityeven observables (m j j , | y j j |, p t, ). However, the interference contribution produces large asymmetric effects in the parity-odd φ j j observable. Constraints are therefore placed on Wilson coefficients using the measured EW Z j j differential cross-section as a function of φ j j .
The measured differential cross-section as a function of φ j j and the corresponding EFT-dependent theoretical prediction are used to define a likelihood function. Statistical correlations amongst the bins of φ j j in the EW Z j j measurement are estimated using a bootstrap procedure (as outlined in Sect. 6) and included in the likelihood function. Each source of systematic uncertainty in the measurement is implemented as a Gaussian-constrained nuisance parameter and is hence treated as fully correlated across bins, but uncorrelated with other uncertainty sources. Uncertainties in the theoretical prediction are also implemented as Gaussian-constrained nuisance parameters. These uncertainties include (i) scale and PDF uncertainties in the Herwig7+Vbfnlo prediction, (ii) an additional shape uncertainty defined by the difference between the Her- operators on the EW Z j j differential cross-sections. The expected contributions from the pure dimension-six term (|M d6 | 2 ) and from the interference between the SM and dimension-six amplitudes (2 Re(M * SM M d6 )) are shown relative to the pure-SM prediction and represented as dotted and dashed lines, respectively. The total contribution to the EW Z j j cross-section is shown as a solid line wig7+Vbfnlo and Powheg+Py8 predictions, and (iii) an uncertainty in the bin-dependent K -factor that arises from finite statistics in the MC samples. The confidence level at each value of Wilson coefficient is calculated using the profile-likelihood test statistic [81], which is assumed to be distributed according to a χ 2 distribution with one degree of freedom following from Wilks' theorem [82]. This allows the 95% confidence intervals to be constructed for each Wilson coefficient. The expected 95% coverage is validated by generating pseudo-experiments, both around the SM hypothesis and at various points in the EFT parameter space.
The expected and observed 95% confidence intervals on the dimension-six operators are shown in Table 4. For each Wilson coefficient, confidence intervals are shown when including or not-including the pure dimension-six contribution in the theoretical prediction. As expected from Fig. 10, the 95% confidence intervals are almost unaffected if the pure dimension-six contributions are excluded from the theoret-ical prediction. The compatibility with the SM hypothesis is found to be poor for one of the operators (Õ H W B ), with a corresponding p-value of 1.6%. The probability that fluctuations around the SM prediction cause this feature when constraining these four Wilson coefficients is investigated using pseudo-experiments. For each pseudo-experiment, the p-value for the compatibility with the SM hypothesis is calculated for each of the four Wilson coefficients. The fraction of pseudo-experiments that produce a p-value lower than 1.6% for any of the Wilson coefficients is found to be 6.2%.
The 95% confidence intervals for the CP-even and CPodd operators can be translated into the HISZ basis [83][84][85] and be compared with previous ATLAS and CMS results.  Table 4 Expected and observed 95% confidence interval for the four Wilson coefficients, using fits to the EW Z j j differential cross-section measured as a function of φ j j . Results are presented when including or excluding the pure dimension-six contributions to the EFT prediction. The p-value quantifying the compatibility with the SM hypothesis is also shown for each Wilson coefficient. The global p-value associated with constraining these four Wilson coefficients is investigated using pseudo-experiments, as outlined in the text tively. These confidence intervals are slightly weaker in sensitivity than the confidence intervals derived using measurements of W + W − production at ATLAS [86], W Z production at CMS [87], and measurements of EW Z j j production at CMS [23]. However, the constraints from those previous measurements were obtained with the pure dimension-six terms included in the theoretical prediction and therefore are more sensitive to the impact of missing higher-dimensional operators in the effective field theory expansion. For example, the constraints obtained from measurements of W W and W Z production are shown to weaken by a factor of ten when the pure dimension-six terms are excluded, due to helicity selection rules that suppress the interference contribution in diboson processes [88,89]. Similarly, the constraints obtained from EW Z j j production at CMS were obtained from a fit to the p t, distribution, which can be dominated by the pure dimension-six terms as shown in Fig. 10. The results presented in this paper therefore have two novel aspects. First, they constitute the strongest limits when pure dimension-six contributions are excluded from the theoretical prediction. Second, the limits are derived from a parity-odd observable, which is sensitive to the interference between the SM and CPodd amplitudes and is therefore a direct test of CP invariance in the weak-boson self-interactions [5].

Conclusion
Differential cross-section measurements for the electroweak production of dijets in association with a Z boson (EW Z j j) are presented for the first time, using proton-proton collision data collected by the ATLAS experiment at a centre-of-mass energy of √ s = 13 TeV and with an integrated luminosity of 139 fb −1 . This process is defined by the t-channel exchange of a weak vector boson and is extremely sensitive to the vector-boson fusion process. Measurements of elec-troweak Z j j production therefore probe the W W Z interaction and provide a fundamental test of the SU(2) L ×U(1) Y gauge symmetry of the Standard Model of particle physics.
The differential cross-sections for EW Z j j production are measured in the Z → + − decay channel ( = e, μ) as a function of four observables: the dijet invariant mass, the rapidity interval spanned by the two jets, the signed azimuthal angle between the two jets, and the transverse momentum of the dilepton pair. The data are corrected for detector inefficiency and resolution using an iterative Bayesian method and are compared to state-of-the-art theoretical predictions from Powheg+Pythia8, Herwig7+Vbfnlo and Sherpa. The data favour the prediction from Herwig7+Vbfnlo . Powheg+Pythia8 predicts too large a cross-section at high values of dijet invariant mass, at large dijet rapidity intervals, and at intermediate values of dilepton transverse momentum. Sherpa predicts too small a cross-section across the measured phase space. Differential cross-section measurements for inclusive Z j j production are also provided in the signal and control regions used to extract the electroweak component.
The detector-corrected measurements are used to search for signatures of anomalous weak-boson self-interactions using the framework of a dimension-six effective field theory. The signed azimuthal angle between the two jets is found to be the most sensitive observable when examining the impact of both the CP-even and CP-odd dimension-six operators. The dimension-six operators are found to be primarily constrained by the contribution to the cross-section from the interference between the SM and dimension-six scattering amplitudes. This makes the results less sensitive to missing higher-order operators in the effective field theory expansion when compared to previous results that search for anomalous weak-boson self-interactions. Furthermore, all limits are derived from a parity-odd observable, which is sensitive to the interference between the SM and CP-odd amplitudes and is therefore a direct test of charge conjugation and parity invariance in the weak-boson self-interactions.

Data Availability Statement
This manuscript has no associated data or the data will not be deposited. [Authors' comment: All ATLAS scientific output is published in journals, and preliminary results are made available in Conference Notes. All are openly available, without restriction on use by external parties beyond copyright law and the standard conditions agreed by CERN. Data associated with journal publications are also made available: tables and data from plots (e.g. cross section values, likelihood profiles, selection efficiencies, cross section limits, ...) are stored in appropriate repositories such as HEPDATA (http:// hepdata.cedar.ac.uk/). ATLAS also strives to make additional material related to the paper available that allows a reinterpretation of the data in the context of new theoretical models. For example, an extended encapsulation of the analysis is often provided for measurements in the framework of RIVET (http://rivet.hepforge.org/). This information is taken from the ATLAS Data Access Policy, which is a public document that can be downloaded from http://opendata.cern.ch/record/413 [opendata.cern.ch].] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 .

A Validation of electroweak extraction methodology
The method used to extract the EW Z j j event yield uses three control regions to constrain the modelling of the strong Z j j background process (see Sect. 5). In this appendix, additional details and validations of the method are presented. First, the EW yields extracted with different strong Z j j predictions is presented. The spread of the extracted yields constitutes the dominant modelling uncertainty in the measured EW Z j j differential cross-sections. Second, EW yields are presented for variations of the electroweak extraction method.
Impact of strong Z j j generator choice Figure 11 shows the EW Z j j event yields extracted using the Sherpa, MG5_NLO+Py8' and MG5+Py8 event generators to predict the strong Z j j background. The results are presented for each bin of the four measured distributions (m j j , | y j j |, p t, and φ j j ). Also shown is the nominal measurement, the central value of which is taken to be the midpoint of the EW Z j j event yields extracted using the three strong Z j j predictions. The EW Z j j yields obtained using Sherpa for the strong Z j j process are typically the largest, while those obtained using MG5+Py8 for the strong Z j j process are typically the lowest. The uncertainty on the nominal measurement due to strong Z j j generator choice is defined as the envelope of the EW Z j j event yields obtained using the three strong Z j j predictions.

Variations of the electroweak extraction method
The degree to which the measured EW Z j j differential crosssections depend on the electroweak extraction method is investigated in this section. First, three variations are applied to the nominal likelihood-based extraction method presented in Sect. 5. These are: 1. Inverted Control Regions: In this variation, CRa and CRc are swapped in Eq. 2. This means that the b i factors are constrained in control regions at high ξ Z and the f (x i ) function is then defined to correct for non-closure when transferring these corrections to low ξ Z . 2. Function choice: In this variation, the choice of f (x i ) in Eq. 2 is changed from a first-order polynomial to a second-order polynomial. Fig. 11 EW Z j j event yields obtained when using the Sherpa, MG5_NLO+Py8' and MG5+Py8 predictions for the strong Z j j process. The extracted EW Z j j yields are presented for each bin of the four observables (m j j , | y j j |, φ j j and p t, ). The nominal EW Z j j measurement is defined as the midpoint of the three generator-specific measurements, with a 'generator choice' systematic defined by the envelope. Experimental systematic uncertainties are not shown (they tend to be small as can be seen in Fig. 6) 3. Unconstrained strong Z j j yield: In this variation, the constraint applied to the integrated strong Z j j event yield (Eq. 3) is removed. This allows the integrated strong Z j j yield to be different between different differential distributions.
In addition, a simpler 'sequential' method is used to extract the EW Z j j yields. In this approach, the EW Z j j yield of bin i in the signal region (SR) is defined by where N data SR,i is the observed event yield in that bin, ν strong,MC SR,i is the strong Z j j background yield predicted by simulation, and ν other,MC SR,i is the event yield from other simulated background processes. The data-driven constraints on the strong Z j j background, r CRa,i and k, are defined by where ν EW,MC CRa,i is the predicted EW Z j j contamination in the control region and is estimated using Herwig7+Vbfnlo .
The index '0' in the definition of k specifies a normalisation region defined by 250 ≤ m j j < 500 GeV. The electroweak signal yield in each bin of the differential distribution is taken to be the midpoint of the envelope of yields obtained using the Sherpa, MG5_NLO+Py8' and MG5+Py8 predictions for the strong Z j j background. The systematic uncertainties on the sequential method are evaluated in the same way as for the nominal method, using the procedures outlined in Sect. 7. Non-closure of the sequential method is evaluated by deriving the constraints on the strong Z j j process using data in CRb and applying them to the strong Z j j prediction in CRc. This non-closure is assigned as an additional systematic uncertainty in the sequential method. Figure 12 shows the EW Z j j event yields measured with the four variations in the electroweak extraction method. The results are found to be in good agreement with the nominal method. In particular, the deviation in the extracted event yields tend to lie within the uncertainty that arises from the strong Z j j generator choice. The agreement between the nominal result and the result obtained with each variation in extraction method is quantified using a χ 2 test for each distribution, with the correlations between measurements determined using the bootstrap method [71]. Good agreement is Fig. 12 Extracted EW Z j j event yields in bins of the four observables measured using five different methods. The nominal measurements with their associated uncertainties are shown as bands, while the measurements based on the four variations to the EW Z j j extraction methodology are shown as points with error bars corresponding to their associated total uncertainty (excluding the experimental systematics, that tend to be significantly smaller) found between the nominal result and each of the results obtained with the different electroweak extraction methods.

B Tabulated differential cross-section measurements
In this section, the measured EW Z j j and inclusive Z j j differential cross-sections are presented in tabular form. The differential cross-sections measured as a function of m j j , | y j j |, p t, , and φ j j are presented in Tables 5-8. The EW Z j j differential cross-sections are measured in the signal region, whereas the inclusive Z j j differential crosssections are measured in the signal region and the three control regions. The fiducial definition for the signal and control regions are defined in Sect. 6. Table 5 Differential cross-section measurements for EW Z j j production and inclusive Z j j production as a function of m j j . The EW Z j j measurements are performed in the signal region. The inclusive Z j j measurements are performed in the signal region and three control regions. The different sources of uncertainty are grouped together for each measurement, with a more granular breakdown of each uncertainty available in HEPDATA EW Z j j SR, m j j cross-section measurements  Table 6 Differential cross-section measurements for EW Z j j production and inclusive Z j j production as a function of | y j j |.