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

Abstract

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 $\sqrt{s}=13\ \hbox {TeV}$ and with an integrated luminosity of $139\ \hbox {fb}^{-1}$. The differential cross-sections are measured in the $Z\rightarrow \ell ^+\ell ^-$ decay channel ($\ell =e,\mu$) 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.

1 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,2,3,4,5,6,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], heavy-vector 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 Zjj production). The EW Zjj 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 Zjj production by ATLAS [19, 20] and CMS [21,22,23] have focused on measuring only an integrated fiducial cross-section in a VBF-enhanced topology. The analysis presented in this article measures differential cross-sections of EW Zjj production in the $Z\rightarrow \ell ^+\ell ^-$ decay channel ($\ell =e,\mu$) and as a function of four observables; the transverse momentum of the dilepton pair ($p_{\textsc {t},\ell \ell }$), the dijet invariant mass ($m_{jj}$), the absolute rapidity¹ separation of the two jets ($|\Delta y_{jj}|$), and the signed azimuthal angle between the two jets ($\Delta \phi _{jj}$). The $\Delta \phi _{jj}$ variable is defined as $\Delta \phi _{jj} = \phi _f -\phi _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 $\sqrt{s}=13\ \hbox {TeV}$ and with an integrated luminosity of $139\ \hbox {fb}^{-1}$.

The EW Zjj 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 Zjj process, the measurements can be used to distinguish between the SM EW Zjj 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 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 $\Delta \phi _{jj}$, for example, is found to be particularly sensitive to anomalous weak-boson self-interactions 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 CP-violating effects in Higgs boson production [24], but has not yet been measured in a final state sensitive to anomalous weak-boson self-interactions.

Fig. 1

Representative Feynman diagrams for EW Zjj production (a, b) and strong Zjj production (c, d). The electroweak Zjj process is defined by the t-channel exchange of a weak boson and at tree level is calculated at $O(\alpha _\text {EW}^4)$ when including the decay of the Z boson. The strong Zjj process has no weak boson exchanged in the t-channel and at tree level is calculated at $O(\alpha _\text {EW}^2\alpha \textsc {s}^2)$ when including the decay of the Z boson

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 Zjj 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 Zjj production are presented in Sect. 8. Differential cross-sections for inclusive Zjj production are also presented in Sect. 8 for the signal and control regions used to extract the electroweak component. The EW Zjj 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.

2 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 $|\eta | < 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 reconstruction up to $|\eta | = 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 $|\eta | < 4.9$. Within the region $|\eta |< 3.2$, electromagnetic calorimetry is provided by barrel and endcap high-granularity lead/liquid-argon (LAr) calorimeters, with an additional thin LAr presampler covering $|\eta | < 1.8$, to correct for energy loss in material upstream of the calorimeters. Hadronic calorimetry is provided by the steel/scintillator-tile calorimeter, segmented into three barrel structures within $|\eta | < 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 $|\eta | < 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 $|\eta | < 2.4$ with resistive-plate chambers in the barrel, and thin-gap chambers in the endcap regions.

Table 1 Summary of generators used for simulation. The details and the corresponding references are provided in the body of the text. In the final column, ‘default’ refers to the default set of tuned parameters provided with the event generator

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.

3 Dataset and Monte Carlo event simulation

The analysis is performed on proton–proton collision data at a centre-of-mass energy of $\sqrt{s}=13\ \hbox {TeV}$. The data were recorded between 2015 and 2018 and correspond to an integrated luminosity of $139\ \hbox {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 Zjj production was simulated using three MC event generators. The default EW Zjj sample was produced with Powheg-Box v1 [29,30,31] using the CT10nlo [32] parton distribution functions (PDF) and is accurate to next-to-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 Zjj production.

The second EW Zjj 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 Zjj production.

The third EW Zjj 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 underlying-event 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 Zjj production.

The dominant background arises from Zjj 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 Zjj background and was simulated using three different MC event generators. Sherpa 2.2.1 was used to produce Z+n-parton 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-next-to-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 Zjj production.

The second strong Zjj 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 Zjj production.

The third strong Zjj sample was also produced with MadGraph5_aMC@NLO, but with the Z+n-parton matrix-elements 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 underlying-event 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 Zjj production.

Production of diboson (VV) 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 OpenLoops 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 ($t\bar{t}$) pair were estimated at NLO accuracy, using the hvq program [54] in Powheg-Box v2. The parton-level 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_{\mathrm{damp}}$ parameter in the Powheg-Box was set to $1.5\,m_{\mathrm{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 (VVV 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 underlying-event 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_{\textsc {t}}$- and $\eta$- 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 Zjj 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 Zjj simulation produced by Powheg+Py8  to match Herwig7+Vbfnlo   at particle level. This is referred to as Herwig7+Vbfnlo’.

4 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_{\textsc {t}}>12~\hbox {GeV}$, whereas the dimuon triggers selected events with leading (subleading) muon candidates having $p_{\textsc {t}}>18\,(8)~\hbox {GeV}$. The transverse momentum thresholds for the lepton candidates were gradually increased during data taking, such that both electron candidates had $p_{\textsc {t}}>24~\hbox {GeV}$ in 2018, whereas the leading muon threshold was increased to 22 GeV in the same running period.

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 Zjj, which is denoted ‘N/A’ in the table. The statistical uncertainty on each prediction is negligible

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 proton–proton interactions are reconstructed using tracking information from the inner detector, with each associated vertex required to have at least two tracks with $p_{\textsc {t}}>0.5~\hbox {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/\sigma _{d_0}| < 3$ and $|z_0 \times \mathrm {sin}\theta | < 0.5~\hbox {mm}$, where $d_0$ is the transverse impact parameter calculated with respect to the measured beam-line position, $\sigma _{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_{\textsc {t}}>25~\hbox {GeV}$ and $|\eta |<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\rightarrow 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/\sigma _{d_0}| < 5$ and $|z_0 \times \mathrm {sin}\theta | < 0.5~\hbox {mm}$. Reconstructed electrons are used in the analysis if they have $p_{\textsc {t}}>25\,\hbox {GeV}$ and $|\eta |<2.47$, but excluding the transition region between the barrel and end-cap calorimeters ($1.37<|\eta |<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 $\gamma$+jet, Z+jet and multijet topologies. Jets are used in the analysis if they have $p_{\textsc {t}}>25\,\hbox {GeV}$ and $|y|<4.4$. As all high-$p_{\textsc {t}}$ electrons pass the above requirements, jets are required to not overlap with a reconstructed electron (i.e. $\Delta R(j,e)>0.2$). Jets with $p_{\textsc {t}}<120\,\hbox {GeV}$ and $|\eta |<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 ($\mathrm {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 Zjj production. A Z-boson candidate is reconstructed by requiring that each event contains exactly two charged leptons ($\ell =e,\mu$) that are opposite in charge and of the same flavour. These leptons are required to be well separated from jets by imposing $\Delta R(\ell ,j)>0.4$. The invariant mass and transverse momentum of the dilepton system is required to fulfil $m_{\ell \ell }\in (81,101)~\hbox {GeV}$ and $p_{\textsc {t},\ell \ell }>20~\hbox {GeV}$. Events are required to contain two or more jets, with the leading and subleading jets satisfying $p_{\textsc {t}}> 85\,\hbox {GeV}$ and $p_{\textsc {t}}> 80\,\hbox {GeV}$, respectively. The dijet system is then constructed from the two leading jets and is required to fulfil $m_{jj}>1\,\hbox {TeV}$ and $|\Delta y_{jj}|>2.0$. The Z boson is required to be centrally produced relative to the dijet system by imposing $\xi _Z<1.0$; the quantity $\xi _Z$ is defined as $\xi _Z =|y_{\ell \ell } - 0.5(y_{j1}+ y_{j2})|\,/\,|\Delta y_{jj}|$, where $y_{\ell \ell }$, $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_{\textsc {t}}^{\mathrm{bal}}<0.15$, where $p_{\textsc {t}}^{\mathrm{bal}}=|\Sigma _i \vec {p}_{\textsc {t},i}|\,/\,\Sigma _i {p}_{\textsc {t},i}$ and the summation includes the dilepton system, the dijet system, and the highest transverse-momentum additional jet reconstructed in the rapidity interval spanned by the dijet system.

Fig. 2

Event yields as a function of $m_{jj}$  (top left), $|\Delta y_{jj}|$  (top right), $p_{\textsc {t},\ell \ell }$  (bottom left) and $\Delta \phi _{jj}$  (bottom right) in data and simulation, measured after the event selection described in Sect. 4. The data are represented as black points and the associated error bar includes only statistical uncertainties. The $m_{jj}$  spectrum is shown starting from 250 GeV, and hence includes more events than the other plots that use the default $m_{jj}> 1000\,\hbox {GeV}$ criterion

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 Zjj event yields predicted by the different event generators. Furthermore, the predicted strong Zjj 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_{jj}$, $|\Delta y_{jj}|$, $p_{\textsc {t},\ell \ell }$, and $\Delta \phi _{jj}$, with Sherpa  used to model the strong Zjj process and Powheg+Py8  used for the EW Zjj process. The level of agreement between data and simulation depends on the kinematic properties of the event, with agreement at large $m_{jj}$  being particularly poor for this configuration of MC simulations.

5 Extraction of electroweak component

Fig. 3

Ratio of Monte Carlo prediction to data for different physics processes and generators for the $m_{jj}$  and $p_{\textsc {t}, \ell \ell }$ 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_{jj}$  distribution extends down to 250 GeV and hence includes a larger phase space than the $p_{\textsc {t},\ell \ell }$ distribution, which requires $m_{jj}>1000\,\hbox {GeV}$. Only statistical uncertainties are shown. The prediction labelled MG5_NLO+Py8’  for the strong Zjj prediction is obtained by a particle-level reweighting of the strong Zjj simulation provided by MG5+Py8. The EW Zjj prediction labelled Herwig7+Vbfnlo’ is also obtained by a particle-level reweighting of the EW Zjj simulation provided by Powheg+Py8  

The poor agreement between data and simulation observed in Fig. 2 implies that the EW Zjj 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 Zjj simulations are used, as shown in Fig. 3 for the $m_{jj}$  and $p_{\textsc {t},\ell \ell }$ distributions. A data-driven method is therefore used to constrain both the shape and normalisation of the strong Zjj background during the extraction of the EW Zjj event yield.

The data are split into four regions by imposing criteria on $\xi _{Z}$   as well as on the multiplicity of jets in the rapidity interval between the leading and subleading jets, $N_{\mathrm{jets}}^{\mathrm{gap}}$. These two variables are chosen because they are almost uncorrelated for both the strong and EW Zjj processes, with calculated correlation coefficients ranging from $-0.04$ to $+0.02$ depending on the event generator and process. Approximately 80% of the EW Zjj events are predicted to fall into the EW-enhanced signal region (SR) defined by $N_{\mathrm{jets}}^{\mathrm{gap}}=0$ and $\xi _{Z}<0.5$. The remaining three regions define EW-suppressed control regions (CR), which can be used to constrain the dominant background from strong Zjj production. These regions are labelled as CRa ($N_{\mathrm{jets}}^{\mathrm{gap}}\ge 1,~\xi _{Z}<0.5$), CRb ($N_{\mathrm{jets}}^{\mathrm{gap}}\ge 1,~\xi _{Z}>0.5$) and CRc ($N_{\mathrm{jets}}^{\mathrm{gap}}=0,~\xi _{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.

Fig. 4

Definition of the signal region (SR) and control regions (CRa, CRb, CRc) used in the extraction of the electroweak component

The EW Zjj event yield is measured in the EW-enhanced SR using a binned maximum-likelihood fit [67, 68]. The log likelihood is defined according to

$$\begin{aligned} \ln {\mathcal {L}} = - \sum _{r,i}\nu _{ri}(\varvec{\theta }) + \sum _{r,i} N_{ri}^{\mathrm{data}}\,\ln {\nu _{ri}(\varvec{\theta })} -\sum _s \frac{\theta _s^2}{2}, \end{aligned}$$

where r is an index corresponding to the region $r\in \left\{ \mathrm {CRa},\mathrm {CRb},\mathrm {CRc},\mathrm {SR}\right\}$, i is the bin of the kinematic observable, $N_{ri}^{\mathrm{data}}$ is the observed event yield and $\nu _{ri}(\varvec{\theta })$ is the prediction that is dependent on the s sources of experimental systematic uncertainty that are each constrained by nuisance parameters $\varvec{\theta } =(\theta _1,\dots , \theta _s$).² The fitted number of events in each region and in each bin of a distribution is given by

$$\begin{aligned} \nu _{ri} = \mu _i\,\nu _{ri}^{\mathrm{EW},\mathrm{MC}} + \nu _{ri}^{\mathrm{strong}} + \nu _{ri}^{\mathrm{other},\mathrm{MC}}, \end{aligned}$$

(1)

where $\mu _i$ is the EW Zjj signal strength of bin i, $\nu _{ri}^{\mathrm{EW},\mathrm{MC}}$ and $\nu _{ri}^{\mathrm{other},\mathrm{MC}}$ are the MC predictions of EW Zjj and contributions from other processes (diboson, $t\bar{t}$ and single top), respectively. The strong Zjj prediction is constrained using the different EW-suppressed control regions according to

$$\begin{aligned} \nu _{\mathrm{CRa},i}^{\mathrm{strong}}= & {} b_{\mathrm{L},i}\,\nu _{\mathrm{CRa},i}^{\mathrm{strong},\mathrm{MC}},\nonumber \\ \qquad \nu _{\mathrm{CRb},i}^{\mathrm{strong}}= & {} b_{\mathrm{H},i}\,\nu _{\mathrm{CRb},i}^{\mathrm{strong},\mathrm{MC}},\nonumber \\ \nu _{\mathrm{SR}i}^{\mathrm{strong}}= & {} b_{\mathrm{L},i}\,f(x_i)\,\nu _{\mathrm{SR},i}^{\mathrm{strong},\mathrm{MC}},\nonumber \\ \qquad \nu _{\mathrm{CRc},i}^{\mathrm{strong}}= & {} b_{\mathrm{H},i}\,f(x_i)\,\nu _{\mathrm{CRc},i}^{\mathrm{strong},\mathrm{MC}}. \end{aligned}$$

(2)

Here, the $b_{\mathrm{L},i}$ and $b_{\mathrm{H},i}$ are sets of bin-dependent factors that apply to the $\xi _{Z}< 0.5$ and $\xi _{Z}>0.5$ regions, respectively. These factors are primarily constrained in CRa and CRb, where they adjust the predicted simulated strong Zjj 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 Zjj yield to account for the extrapolation from CRa ($N_{\mathrm{jets}}^{\mathrm{gap}}\ge 1$) to the SR ($N_{\mathrm{jets}}^{\mathrm{gap}}=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 $\mu _i$, the two parameters of the function $f(x_i)$, and the $b_{\mathrm{L},i}$ and $b_{\mathrm{H},i}$ corrections to the strong Zjj process. In total, this amounts to $3\,N_{\mathrm{bins}}+2$ parameters that are constrained using $4\,N_{\mathrm{bins}}$ measurements in data, where $N_{\mathrm{bins}}$ is the number of bins measured for a specific observable ($m_{jj}$, $|\Delta y_{jj}|$, $p_{\textsc {t}, \ell \ell }$  and $\Delta \phi _{jj}$).

The pre-fit and post-fit agreement between data and simulation is shown in Fig. 5 as a function of $m_{jj}$  in the signal and control regions. Two separate fits are shown, one using the Sherpa  strong Zjj 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_{jj}$, but produce very good agreement with the data following the fitting procedure. The overall scaling factor applied to the strong Zjj prediction from MG5_NLO+Py8’   in the signal region is 0.93 at low $m_{jj}$  rising to 2.2 at high $m_{jj}$. For Sherpa, the corresponding scaling factors are 0.86 at low $m_{jj}$   and 0.26 at high $m_{jj}$. The pre-fit systematic uncertainties shown on the plots are derived as outlined in Sect. 7.

Fig. 5

Comparison between data and prediction before (left) and after (right) the fit using strong Zjj estimates based on Sherpa  (top) and MG5_NLO+Py8’  (bottom) in bins of $m_{jj}$  in the different control and signal regions. The MG5_NLO+Py8’  prediction is obtained by a particle-level reweighting of the strong Zjj simulation provided by MG5+Py8. The $m_{jj}$  bin edges are defined by $(1.0,1.5,2.25,3.0,4.5,7.5)\ \hbox {TeV}$

Since there is no a priori reason to prefer any strong Zjj generator over another, the EW Zjj component is extracted three times, once using the Sherpa  strong Zjj prediction, once using the MG5_NLO+Py8’   strong Zjj prediction, and once using the MG5+Py8   strong Zjj 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 Zjj event generators. The envelope itself is used to define a systematic uncertainty as outlined in Sect. 7.

The constraints on the strong Zjj simulation in Eq. 2 are evaluated independently for each of the measured differential distributions ($m_{jj}$, $|\Delta y_{jj}|$, $p_{\textsc {t},\ell \ell }$   and $\Delta \phi _{jj}$). This results in slightly different total EW Zjj and strong Zjj 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 Zjj yield is obtained for each distribution, i.e.

$$\begin{aligned} \sum _i \nu _{\mathrm{SR},i}^{\mathrm{strong}} = \hat{\nu }^{\mathrm{strong}}_{\mathrm{SR},m_{jj}}, \end{aligned}$$

(3)

where $\hat{\nu }^{\mathrm{strong}}_{\mathrm{SR},m_{jj}}$ is the event yield obtained by integrating the constrained strong Zjj template for the $m_{jj}$  distribution in the SR .

Table 3 Particle-level definition of the measurement. $\Delta R_{\mathrm{min}}(\ell _1,j)$ denotes the minimum $\Delta R$ distance between the highest transverse-momentum lepton ($\ell _1$) and any of the jets in the event. $\Delta R_{\mathrm{min}}(\ell _2,j)$ is similarly defined

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 Zjj simulation as defined in Eq. 2, such that the $b_i$ factors are constrained in CRs at high $\xi _Z$ and the $f(x_i)$ function is then defined to correct for non-closure when transferring these corrections to low $\xi _Z$. Second, the constraint on the strong Zjj 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 second-order polynomial. Third, the constraint applied to the integrated strong Zjj event yield (Eq. 3) is removed. Finally, a simpler ‘sequential’ method is used to extract the EW Zjj event yields. In this approach, the data-driven correction to the strong Zjj is derived in CRa (assuming the SM prediction for the electroweak process in this region) and directly applied to the strong Zjj simulation in the SR. A transfer factor to account for mismodelling between the SR and CRa is evaluated at low $m_{jj}$   ($250\le m_{jj}<500~\hbox {GeV}$). Non-closure of the sequential method is evaluated in CRc using corrections to the strong Zjj 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.

6 Correction for detector effects

Particle-level differential cross-sections are produced by correcting the inclusive Zjj and EW Zjj event yields in each bin for the effects of detector inefficiency and resolution. The EW Zjj event yields are extracted in the signal region using the method outlined in the previous section. The inclusive Zjj 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 ($t\bar{t}$, single-top, $VV\not \rightarrow Zjj$, and W+jets production). For both inclusive and EW Zjj production, the event yields in the $e^+e^-$ and $\mu ^+\mu ^-$ 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\tau > 10~\hbox {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 $\Delta R < 0.1$. Leptons are required to have $p_{\textsc {t}}>25~\hbox {GeV}$ and have the same acceptance requirement as used at the analysis level, i.e. muons satisfy $|\eta |<2.4$ and electrons satisfy $|\eta |<2.47$ (but exclude the region $1.37<|\eta |<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_{\textsc {t}}>25\,\hbox {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 Zjj differential cross-sections are measured in the SR, whereas inclusive Zjj 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 detector-level selection. For the EW Zjj differential cross-section measurements, the Powheg+Py8   EW Zjj simulation is used to define the corrections and the response matrices. For the inclusive Zjj differential cross-section measurements, all sources of Zjj production are part of the measurement and the unfolding is carried out using the cross-section weighted sum of the Powheg+Py8   EW Zjj simulation, the Sherpa  strong Zjj 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 Zjj measurements, the electroweak extraction is repeated for each pseudo-experiment after fluctuating the event yields, in each bin of the signal and control regions, using a Poisson distribution. For the inclusive Zjj measurements, the background-subtracted 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.

7 Systematic uncertainties

7.1 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 Zjj process, as well as the MC simulations used to unfold the EW Zjj and inclusive Zjj 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 four-momentum 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_{jj}$   or $p_{\textsc {t},\ell \ell }$, but more than 10% for $m_{jj}>4\,\hbox {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_{\mathrm{jets}}^{\mathrm{gap}}=0$ and $N_{\mathrm{jets}}^{\mathrm{gap}}\ge 1$ regions.

7.2 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 Zjj process are the dominant uncertainties in the extraction of the electroweak signal yield. Three sources of uncertainty in the strong Zjj modelling are investigated, arising from (i) the choice of event generator, (ii) the renormalisation and factorisation scale dependence in the strong Zjj 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 Zjj 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 Zjj templates for Sherpa  in which the renormalisation ($\mu _{{\textsc {R}}}$) and factorisation ($\mu _{{\textsc {F}}}$) scales have been varied independently by factors of 0.5 and 2.0. Six variations are considered for each generator corresponding to $(\mu _{{\textsc {R}}}, \mu _{{\textsc {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 Zjj 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 Zjj sample to reproduce the variations of the NNPDF3.0nnlo PDF set (including the associated $\alpha \textsc {s}$variations) 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 Zjj production, the choice of event generator has the largest impact on the extracted electroweak yields.

Theoretical uncertainties associated with the modelling of the EW Zjj 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 Zjj process is determined by repeating the analysis after reweighting the default Powheg-Box EW Zjj simulation such that it matches the prediction of the Herwig7+Vbfnlo   EW Zjj 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 Zjj simulation to extract the electroweak event yield and the results are found to be consistent and within the assigned uncertainty due to electroweak Zjj 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 Zjj 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 Zjj process and the strong Zjj 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 Zjj 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 Zjj process.

Fig. 6

Fractional uncertainty in the inclusive Zjj measurement (top) and the EW Zjj measurement (bottom) as a function of $m_{jj}$   (left) and $p_{\textsc {t}, \ell \ell }$    (right). Uncertainty sources are grouped in categories that are added in quadrature (denoted $\oplus$) to give the total uncertainty. The ‘EW Zjj model’ component includes the uncertainty on the EW Zjj prediction and the impact of interference between the strong Zjj and EW Zjj processes. The ‘strong Zjj model’ uncertainty is dominated by the choice of generator used for the strong Zjj prediction, but also includes the impact of renormalisation/factorisation scale variations and PDF set variations

7.3 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 Zjj differential cross-section measurements, the Sherpa  EW Zjj simulation is used in place of the Powheg+Py8  EW Zjj simulation. For the inclusive Zjj differential cross-section measurements, the MG5+Py8  strong Zjj simulation is used in place of the Sherpa  strong Zjj 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 Zjj 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.

7.4 Summary of systematic uncertainties

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

8 Results

The differential cross-sections for EW Zjj production as a function of $m_{jj}$, $|\Delta y_{jj}|$, $p_{\textsc {t}, \ell \ell }$, and $\Delta \phi _{jj}$   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 Herwig7+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 Zjj 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 PDF-related uncertainties on the EW Zjj predictions are found to be much smaller than the impact of scale uncertainties.

Fig. 7

Differential cross-sections for EW Zjj production as a function of $m_{jj}$  (top left), $|\Delta y_{jj}|$  (top right), $p_{\textsc {t}, \ell \ell }$ (bottom left) and $\Delta \phi _{jj}$  (bottom right). 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

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 to overestimate the EW Zjj cross-section at high $m_{jj}$, high $|\Delta y_{jj}|$, and intermediate $p_{\textsc {t}, \ell \ell }$. 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 VVjj 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 Zjj 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 Zjj 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.

Fig. 8

Differential cross-sections measured in the SR for inclusive Zjj production as a function of $m_{jj}$  (top left), $|\Delta y_{jj}|$   (top right), $p_{\textsc {t}, \ell \ell }$ (bottom left) and $\Delta \phi _{jj}$  (bottom right). 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 Zjj 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

A fiducial cross-section for EW Zjj production is calculated, by integrating the differential cross-section as a function of $m_{jj}$, and found to be

$$\begin{aligned} \sigma _{\mathrm{EW}} = 37.4 \pm 3.5 \,(\text {stat}) \pm 5.5\,(\text {syst})\, \, \text {fb}. \end{aligned}$$

This is in excellent agreement with the theoretical prediction from Herwig7+Vbfnlo , which is $39.5 \pm 3.4\,(\text {scale}) \pm 1.2 \,(\text {PDF})~\hbox {fb}$.

Fig. 9

Differential cross-sections measured in CRa for inclusive Zjj production as a function of $m_{jj}$  (left) and $p_{\textsc {t}, \ell \ell }$ (right), where CRa is defined by $N_{\mathrm{jets}}^{\mathrm{gap}}\ge 1$ and $\xi _{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 Zjj 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

Differential cross-sections for inclusive Zjj production as a function of $m_{jj}$, $|\Delta y_{jj}|$, $p_{\textsc {t},\ell \ell }$ and $\Delta \phi _{jj}$  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 Zjj background presumably will become available. The differential cross-sections for inclusive Zjj production measured in the SR as a function of $m_{jj}$, $|\Delta y_{jj}|$, $p_{\textsc {t}, \ell \ell }$, and $\Delta \phi _{jj}$  are shown in Fig. 8. The differential cross-sections measured in CRa for inclusive Zjj production as a function of $m_{jj}$  and $p_{\textsc {t}, \ell \ell }$ are shown in Fig. 9. The data are compared with the strong Zjj predictions provided by Sherpaand MG5_NLO+Py8, augmented with the EW Zjj contribution predicted by Herwig7+Vbfnlo   and the VZ contribution predicted by Sherpa. The effects of scale uncertainties on the strong Zjj 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 Zjj predictions are estimated using the variations of the NNPDF PDF sets. The total uncertainty on the strong Zjj 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 Zjj production.

The unfolded differential cross-sections for EW Zjj production and inclusive Zjj 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].

9 Constraints on anomalous weak-boson self-interactions

In this section, the measured EW Zjj differential cross-sections 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

$$\begin{aligned} {\mathcal {L}}_{\mathrm{eff}} = {\mathcal {L}}_{\mathrm{SM}} + \sum _i \frac{c_i}{\Lambda ^2} {\mathcal {O}}_i , \end{aligned}$$

(4)

where ${\mathcal {L}}_{\mathrm{SM}}$ is the SM Lagrangian, the ${\mathcal {O}}_i$ are dimension-six operators in the Warsaw basis [78], and the ${c_i}/\Lambda ^2$ are Wilson coefficients that describe the strength of the anomalous interactions induced by those operators. Constraints are placed on two CP-even operators (${\mathcal {O}}_W$, ${\mathcal {O}}_{HWB}$) and two CP-odd operators ($\tilde{{\mathcal {O}}}_W$, $\tilde{{\mathcal {O}}}_{HWB}$), which are known to produce anomalous WWZ interactions.

Theoretical predictions are constructed for the EW Zjj process using the effective Lagrangian in Eq. 4. The amplitude for the EW Zjj process is split into a SM part, ${\mathcal {M}}_{\mathrm{SM}}$, and a dimension-six part, ${\mathcal {M}}_{\mathrm{d}6}$, which contains the anomalous interactions. The differential cross-section or squared amplitude then has three contributions

$$\begin{aligned} |{\mathcal {M}}|^2 = |{\mathcal {M}}_{\mathrm{SM}}|^2 + 2\,\mathrm {Re}({\mathcal {M}}^*_{\mathrm{SM}} {\mathcal {M}}^{\,}_{\mathrm{d6}})+|{\mathcal {M}}_{\mathrm{d6}}|^2, \end{aligned}$$

(5)

namely a pure SM term $|{\mathcal {M}}_{\mathrm{SM}}|^2$, a pure dimension-six term $|{\mathcal {M}}_{\mathrm{d6}}|^2$, and a term that contains the interference between the SM and dimension-six amplitudes, $2\,\mathrm {Re}({\mathcal {M}}^*_{\mathrm{SM}} {\mathcal {M}}^{\,}_{\mathrm{d6}})$. The constraints on the dimension-six operators presented in this section are derived both with and without the pure dimension-six terms included in the theoretical prediction. This tests whether the results are robust against missing dimension-eight operators in the EFT expansion.

The pure-SM contribution to the EW Zjj differential cross-sections in Eq. 5 is taken to be the prediction from Herwig7+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 higher-order QCD corrections, the interference and pure dimension-six contributions are scaled using a bin-dependent K-factor, which is defined by the ratio of pure-SM EW Zjj differential cross-sections predicted by Herwig7+Vbfnlo   and MadGraph5$+$ Pythia8  in each bin.

The impact of the interference and pure dimension-six contributions to the EW Zjj differential cross-sections is shown relative to the pure SM contribution in Fig. 10. The Wilson coefficients were chosen to be ${c_W}/\Lambda ^2 = 0.2\ \hbox {TeV}^{-2}$, ${\tilde{c}_W}/\Lambda ^2 = 0.2\ \hbox {TeV}^{-2}$, ${c_{HWB}}/\Lambda ^2 = 1.8\ \hbox {TeV}^{-2}$ and ${\tilde{c}_{HWB}}/\Lambda ^2 = 1.8\ \hbox {TeV}^{-2}$. For the CP-even ${\mathcal {O}}_W$ operator, the high-$p_{\textsc {t},\ell \ell }$ region is particularly sensitive to the anomalous interactions, a feature that was seen in previous studies for EW Vjj production [23, 80]. The pure dimension-six contributions to the cross-section dominate in this region. The $\Delta \phi _{jj}$  observable is also found to be very sensitive to the anomalous interactions induced by the ${\mathcal {O}}_W$ operator, but in this observable the interference contribution dominates. For the CP-even ${\mathcal {O}}_{HWB}$ operator, the interference contribution dominates in all distributions, with the $\Delta \phi _{jj}$   observable showing the largest kinematic dependences. For the CP-odd operators, the interference contribution is zero in the parity-even observables ($m_{jj}$, $|\Delta y_{jj}|$, $p_{\textsc {t}, \ell \ell }$). However, the interference contribution produces large asymmetric effects in the parity-odd $\Delta \phi _{jj}$  observable. Constraints are therefore placed on Wilson coefficients using the measured EW Zjj differential cross-section as a function of $\Delta \phi _{jj}$.

Fig. 10

Impact of the ${\mathcal {O}}_W$, $\tilde{{\mathcal {O}}}_W$, ${\mathcal {O}}_{HWB}$ and $\tilde{{\mathcal {O}}}_{HWB}$ operators on the EW Zjj differential cross-sections. The expected contributions from the pure dimension-six term ($|{\mathcal {M}}_{\mathrm{d6}}|^2$) and from the interference between the SM and dimension-six amplitudes ($2\,\mathrm {Re}({\mathcal {M}}^*_{\mathrm{SM}} {\mathcal {M}}^{\,}_{\mathrm{d6}})$) are shown relative to the pure-SM prediction and represented as dotted and dashed lines, respectively. The total contribution to the EW Zjj cross-section is shown as a solid line

The measured differential cross-section as a function of $\Delta \phi _{jj}$  and the corresponding EFT-dependent theoretical prediction are used to define a likelihood function. Statistical correlations amongst the bins of $\Delta \phi _{jj}$  in the EW Zjj 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 Herwig7+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 $\chi ^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.

Table 4 Expected and observed 95% confidence interval for the four Wilson coefficients, using fits to the EW Zjj differential cross-section measured as a function of $\Delta \phi _{jj}$. 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

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 theoretical prediction. The compatibility with the SM hypothesis is found to be poor for one of the operators ($\tilde{{\mathcal {O}}}_{HWB}$), 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 CP-odd operators can be translated into the HISZ basis [83,84,85] and be compared with previous ATLAS and CMS results. The observed and expected 95% confidence intervals for the $c_{WWW}/\Lambda ^2$ Wilson coefficient are [–2.7, 5.8] $\hbox {TeV}^{-2}$ and [–4.4, 4.1] $\hbox {TeV}^{-2}$, respectively. The observed and expected 95% confidence intervals for the $\tilde{c}_{WWW}/\Lambda ^2$ Wilson coefficient are [–1.6, 2.0] $\hbox {TeV}^{-2}$ and [–1.7, 1.7] $\hbox {TeV}^{-2}$ respectively. These confidence intervals are slightly weaker in sensitivity than the confidence intervals derived using measurements of $W^+W^-$ production at ATLAS  [86], WZ production at CMS [87], and measurements of EW Zjj 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 WW and WZ 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 Zjj production at CMS were obtained from a fit to the $p_{\textsc {t}, \ell \ell }$ 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 CP-odd amplitudes and is therefore a direct test of CP invariance in the weak-boson self-interactions [5].

10 Conclusion

Differential cross-section measurements for the electroweak production of dijets in association with a Z boson (EW Zjj) are presented for the first time, using proton–proton collision data collected by the ATLAS experiment at a centre-of-mass energy of $\sqrt{s}=13\ \hbox {TeV}$ and with an integrated luminosity of $139~\hbox {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 electroweak Zjj production therefore probe the WWZ interaction and provide a fundamental test of the SU(2)$_{\mathrm{L}} \times$U(1)$_{\mathrm{Y}}$ gauge symmetry of the Standard Model of particle physics.

The differential cross-sections for EW Zjj production are measured in the $Z\rightarrow \ell ^+\ell ^-$ decay channel ($\ell =e,\mu$) 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 Zjj 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.

Appendices

A Validation of electroweak extraction methodology

The method used to extract the EW Zjj event yield uses three control regions to constrain the modelling of the strong Zjj 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 Zjj predictions is presented. The spread of the extracted yields constitutes the dominant modelling uncertainty in the measured EW Zjj differential cross-sections. Second, EW yields are presented for variations of the electroweak extraction method.

1.1 Impact of strong Zjj generator choice

Figure 11 shows the EW Zjj event yields extracted using the Sherpa, MG5_NLO+Py8’  and MG5+Py8  event generators to predict the strong Zjj background. The results are presented for each bin of the four measured distributions ($m_{jj}$, $|\Delta y_{jj}|$, $p_{\textsc {t}, \ell \ell }$ and $\Delta \phi _{jj}$). Also shown is the nominal measurement, the central value of which is taken to be the midpoint of the EW Zjj event yields extracted using the three strong Zjj predictions. The EW Zjj yields obtained using Sherpa  for the strong Zjj process are typically the largest, while those obtained using MG5+Py8  for the strong Zjj process are typically the lowest. The uncertainty on the nominal measurement due to strong Zjj generator choice is defined as the envelope of the EW Zjj event yields obtained using the three strong Zjj predictions.

Fig. 11

EW Zjj event yields obtained when using the Sherpa, MG5_NLO+Py8’  and MG5+Py8  predictions for the strong Zjj process. The extracted EW Zjj yields are presented for each bin of the four observables ($m_{jj}$, $|\Delta y_{jj}|$, $\Delta \phi _{jj}$  and $p_{\textsc {t}, \ell \ell }$). The nominal EW Zjj 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)

1.2 Variations of the electroweak extraction method

The degree to which the measured EW Zjj differential cross-sections 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 $\xi _Z$ and the $f(x_i)$ function is then defined to correct for non-closure when transferring these corrections to low $\xi _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.

• 3. Unconstrained strong Zjj yield: In this variation, the constraint applied to the integrated strong Zjj event yield (Eq. 3) is removed. This allows the integrated strong Zjj yield to be different between different differential distributions.

In addition, a simpler ‘sequential’ method is used to extract the EW Zjj yields. In this approach, the EW Zjj yield of bin i in the signal region (SR) is defined by

$$\begin{aligned} \nu _{\mathrm{SR},i}^{\mathrm{EW}} = N_{\mathrm{SR},i}^{\mathrm{data}} - k \, r^{\,}_{\mathrm{CRa},i} \, \nu _{SR,i}^{\mathrm{strong},\mathrm{MC}} - \nu _{\mathrm{SR},i}^{\mathrm{other},\mathrm{MC}}, \end{aligned}$$

(6)

where $N_{\mathrm{SR},i}^{\mathrm{data}}$ is the observed event yield in that bin, $\nu _{\mathrm{SR},i}^{\mathrm{strong},\mathrm{MC}}$ is the strong Zjj background yield predicted by simulation, and $\nu _{\mathrm{SR},i}^{\mathrm{other},\mathrm{MC}}$ is the event yield from other simulated background processes. The data-driven constraints on the strong Zjj background, $r^{\,}_{\mathrm{CRa},i}$ and k, are defined by

$$\begin{aligned}&r^{\,}_{\mathrm{CRa}, i} = \left( N_{\mathrm{CRa},i}^{\mathrm{data}} - \nu _{\mathrm{CRa},i}^{\mathrm{EW},\mathrm{MC}} - \nu _{\mathrm{CRa},i}^{\mathrm{other}, \mathrm{MC}}\right) /\nu _{\mathrm{CRa},i}^{\mathrm{strong},\mathrm{MC}}\\&\qquad \text {and} \quad k= r^{\,}_{\mathrm{SR},0}/r^{\,}_{\mathrm{CRa},0}, \end{aligned}$$

where $\nu _{\mathrm{CRa},i}^{\mathrm{EW},\mathrm{MC}}$ is the predicted EW Zjj 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\le m_{jj}<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 Zjj 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 Zjj process using data in CRb and applying them to the strong Zjj prediction in CRc. This non-closure is assigned as an additional systematic uncertainty in the sequential method.

Fig. 12

Extracted EW Zjj 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 Zjj 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)

Figure 12 shows the EW Zjj 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 Zjj generator choice. The agreement between the nominal result and the result obtained with each variation in extraction method is quantified using a $\chi ^2$ test for each distribution, with the correlations between measurements determined using the bootstrap method [71]. Good agreement is 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 Zjj and inclusive Zjj differential cross-sections are presented in tabular form. The differential cross-sections measured as a function of $m_{jj}$, $|\Delta y_{jj}|$, $p_{\textsc {t}, \ell \ell }$, and $\Delta \phi _{jj}$  are presented in Tables 5–8. The EW Zjj differential cross-sections are measured in the signal region, whereas the inclusive Zjj differential cross-sections 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 Zjj production and inclusive Zjj production as a function of $m_{jj}$. The EW Zjj measurements are performed in the signal region. The inclusive Zjj 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

Table 6 Differential cross-section measurements for EW Zjj production and inclusive Zjj production as a function of $|\Delta y_{jj}|$. The EW Zjj measurements are performed in the signal region. The inclusive Zjj 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

Table 7 Differential cross-section measurements for EW Zjj production and inclusive Zjj production as a function of $p_{\textsc {t},\ell \ell }$. The EW Zjj measurements are performed in the signal region. The inclusive Zjj 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

Table 8 Differential cross-section measurements for EW Zjj production and inclusive Zjj production as a function of $\Delta \phi _{jj}$. The EW Zjj measurements are performed in the signal region. The inclusive Zjj 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