1 Introduction

Leptoquarks (LQs) are features of a number of extensions of the Standard Model (SM) [1,2,3,4,5,6,7,8] and may provide an explanation for the similarities between the quark and lepton sectors in the SM. They also appear in models addressing some of the recent b-flavour anomalies [9,10,11]. They are colour-triplet bosons with fractional electric charge and possess non-zero baryon and lepton number [12]. Scalar and vector LQs have been proposed and are expected to decay directly into lepton–quark pairs. The lepton can be either electrically charged or neutral.

A single Yukawa coupling, \(\lambda _{\mathrm {LQ} \rightarrow \ell q}\), determines the coupling strength between scalar LQs and the lepton–quark pair [13]. Two additional coupling constants due to magnetic moment and electric quadrupole moment interactions are needed for vector LQs [14].

LQs can be produced both singly and in pairs in proton–proton (pp) interactions; diagrams showing representative single- and pair-production processes are shown in Fig. 1. The single-LQ production cross-section depends on \(\lambda _{\mathrm {LQ} \rightarrow \ell q}\) whereas the pair-production cross-section is largely insensitive to this coupling. For pp interactions with a centre-of-mass energy \(\sqrt{s}\) = 13 TeV, gluon fusion represents the dominant pair-production mechanism for LQ masses below around 1 TeV. The contribution of the \(q\bar{q}\)-annihilation process rises with LQ mass. Only scalar LQ production is considered in this paper because this is less model dependent than vector LQ production. The production of vector LQs depends on additional parameters and a full interpretation in such models is beyond the scope of this analysis. The results obtained here can, however, be regarded as conservative estimates of limits on vector LQ production, since the production cross-section for vector LQs is typically much larger than for scalar LQs, while the kinematic properties used to search for their signature are very similar for both spin hypotheses [8].

Fig. 1
figure 1

Diagrams for representative a single and b pair production of LQs

The benchmark signal model used in this paper is the minimal Buchmüller–Rückl–Wyler (BRW) model [15]. This model imposes a number of constraints on the properties of the LQs. Couplings are purely chiral and LQs belong to three families, corresponding to the three SM generations, such that only leptons and quarks within a given generation can interact. The requirement of same-generation interactions excludes flavour-changing neutral currents (FCNC) [16]. The branching ratio (BR) of a LQ decay into different states is taken as a free parameter. In this paper, \(\beta \) denotes the BR for a LQ decay into a charged lepton and quark, \(\mathrm {LQ}\rightarrow \ell ^{\pm }q\). The BR for a LQ decay into a neutrino and quark is \((1 - \beta )\).

This paper presents a search for LQs pair-produced in pp collisions at \(\sqrt{s}\) = 13 TeV, performed by the ATLAS experiment at the Large Hadron Collider (LHC). In previous searches for pair-produced LQs made by ATLAS [17,18,19,20,21], most recently with 3.2 fb\(^{-1}\) of data collected at \(\sqrt{s}\) = 13 TeV, the existence of scalar LQs with masses below 1100 (900) GeV for first-generation LQs at \(\beta \) = 1 (0.5) and, for second-generation LQs, with masses below 1050 (830) GeV at \(\beta \) = 1 (0.5) was excluded [21]. Searches have also been made by the CMS Collaboration [22,23,24,25,26,27,28] using the same model assumptions as in this work. That experiment found that scalar LQs with masses below 1435 (1270) GeV for first-generation LQs at \(\beta \) = 1 (0.5) and below 1530 (1285) GeV for second-generation LQs at \(\beta \) = 1 (0.5) are excluded at 95% confidence level (CL) using a data sample of 35.9 fb\(^{-1}\) collected at \(\sqrt{s}\) = 13 TeV [27, 28]. An overview of limits on LQ production and masses can be found in Ref. [29].

The decay of pair-produced LQs can lead to final states containing a pair of charged leptons and jets, or one charged lepton, a neutrino and jets. For down-type leptoquarks with a charge of \(\frac{1}{3}e\), assumed here as a benchmark, this can be written more explicitly as \(\mathrm {LQ}\overline{\mathrm {LQ}}\longrightarrow \ell ^{-}q\ell ^{+}\bar{q}\) for the dilepton channel, and \(\mathrm {LQ}\overline{\mathrm {LQ}}\longrightarrow \ell ^{-}q\bar{\nu }\bar{q'}\) or \(\mathrm {LQ}\overline{\mathrm {LQ}}\longrightarrow \nu q\ell ^{+}\bar{q'}\) for the lepton-neutrino channel, respectively. The four search channels in this paper correspond to the above decays. The eejj (\(\mu \mu \)jj) channel comprises exactly two electrons (muons) and at least two jets, and the e \(\nu \)jj (\(\mu \nu \)jj) channel requires exactly one electron (muon), missing momentum in the transverse plane and at least two jets. The electron and muon channels are treated separately as they correspond to LQs of different generations. The dilepton and lepton–neutrino channels are combined in order to increase the sensitivity to the parameter of the BRW model, \(\beta \). While the main target of this paper are first- and second-generation LQs, there are also dedicated searches looking for 3rd generation LQ production [20, 30,31,32,33,34].

The main background processes in this search are \(Z/\gamma ^*{+}\)jets, \(W{+}\)jets and \(t\bar{t}\) production. The shapes of these background distributions are estimated from Monte Carlo (MC) simulations while their normalisation is obtained from a simultaneous background-only profile likelihood fit to three data control regions (CRs). The CRs are defined with specific selections such that they are orthogonal to the signal regions (SRs), free from any LQ signal of interest, and enriched in their relevant background process. For the particular benchmark model chosen, further signal and background discrimination could be achieved by exploiting flavour tagging in order to reject events with b-hadrons in the search region. This is not used, in order to maintain sensitivity to a \(\mathrm {LQ}\) coupling to a b-quark. The shapes and normalisations of other (subdominant) backgrounds are estimated from MC simulation. Background contributions from misidentified leptons are also sub-dominant, and are estimated from data when their contribution is significant.

Several variables that provide separation between signal and background processes are combined into a single discriminant using a boosted decision tree (BDT) method [35]. For each LQ-mass hypothesis and for each lepton flavour, two BDTs are used: one for the dilepton channel and one for the lepton–neutrino channel. The BDT output score is binned and the signal is accumulated at high score values. A single bin at a high score defines the SR for the respective mass hypothesis. The bin range is chosen such that a sensitivity estimator based on the expected number of signal and background events is maximised.

A profile likelihood fit is performed simultaneously to both the dilepton and lepton-neutrino channel for each lepton flavour to extract a limit on the signal yield. This is done for different assumptions of LQ mass and branching ratio \(\beta \), resulting in exclusion bounds in terms of these parameters.

In addition to the search results, this paper presents a set of particle-level fiducial and differential cross-section measurements in six dilepton–dijet regions which are related to the analysis CRs. The measurements are made for the dominant process in each region (\(Z{+}\)jets or \(t\bar{t}\)) exclusively. The motivation for making these measurements is as follows. Search CRs are typically in kinematic regions which are not directly covered by dedicated measurements of SM processes. For instance, the ATLAS measurement of \(Z+\text {jets}\) production described in Ref. [36] has much looser selections on the transverse momentum (\(p_{\mathrm {T}}\)) of leptons and jets than the region that includes Z bosons for this search. Event generators are tuned and validated to ensure that they agree with the available SM data. However, in some cases, although good agreement is seen for the regions covered by dedicated SM measurements, significant disagreements remain possible in more extreme regions where no measurements are available to validate the generator predictions. Indeed, although no disagreement was seen when validating generated \(Z{+}\)jets samples against various SM measurements [37], significant disagreement was observed in the \(Z{+}\)jets CR for this search. Extracting particle-level cross-sections in these CRs and related more extreme regions of phase space can therefore provide complementary information to validate the performance of event generators. Such measurements may lead to improved background modelling in future searches with final states using similar CRs (see for example Ref. [38]).

The paper is organised as follows. Section 2 gives an overview of the ATLAS detector, while Sect. 3 summarises information on the data set and simulations used. The physics objects used in the analysis are defined in Sect. 4. The background estimation and CRs used in the process are described in Sect. 5, and Sect. 6 provides details on the cross-section extraction. Section 7 gives an overview over the multivariate analysis employed for the search and the SRs defined with it. The systematic uncertainties are summarised in Sect. 8, and the statistical procedure employed for the search is described in Sect. 9. Results from both the measurement and the search are given in Sect. 10. Section 11 concludes the paper.

2 The ATLAS detector

The ATLAS experiment [39] is a multipurpose detector with a forward–backward symmetric cylindrical geometry and nearly 4\(\pi \) coverage in solid angle.Footnote 1 The three major subcomponents are the tracking detector, the calorimeter and the muon spectrometer. Charged-particle tracks and vertices are reconstructed by the inner detector (ID) tracking system, comprising silicon pixel and microstrip detectors covering the pseudorapidity range \(|\eta |\) < 2.5, and a transition radiation straw-tube tracker that covers \(|\eta |\) < 2.0 and provides electron identification. The ID is immersed in a homogeneous 2 T magnetic field provided by a solenoid. Electron, photon, and jet energies are measured with sampling calorimeters. The calorimeter system covers a pseudorapidity range of \(|\eta |\) < 4.9. Within the region \(|\eta |\) < 3.2, barrel and endcap high-granularity lead/liquid argon (LAr) electromagnetic calorimeters are deployed, 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 a steel/scintillator-tile calorimeter, segmented into three barrel structures within \(|\eta |\) < 1.7, and two copper/LAr hadronic endcap calorimeters. The forward region (3.1 < \(|\eta |\) < 4.9) is instrumented with a LAr calorimeter with copper (electromagnetic) and tungsten (hadronic) absorbers. Surrounding the calorimeters is a muon spectrometer (MS) with air-core toroid magnets to provide precise muon identification and momentum measurements, consisting of a system of precision tracking chambers providing coverage over \(|\eta |\) < 2.7, and detectors with triggering capabilities over \(|\eta |\) < 2.4. A two-level trigger system [40], the first level using custom hardware and followed by a software-based level, is used to reduce the event rate to a maximum of around 1 kHz for offline storage.

3 Data and Monte Carlo samples

3.1 Data sample

This analysis is based on proton–proton collision data at a centre-of-mass energy of \(\sqrt{s}\) = 13 TeV, collected at the LHC during 2015 and 2016. After imposing requirements based on beam and detector conditions and data quality, the data sample corresponds to an integrated luminosity of \(36.1\,{\mathrm{fb}}^{-1}\). A set of triggers is used to select events [40,41,42]. For the dielectron channel, a two-electron trigger with a transverse energy (\(E_{\mathrm {T}}\)) threshold of 17 GeV for each electron was used for the entire data set. In the electron–neutrino channel, events were selected by either of two single-electron triggers, as described below. For the 2015 data set, one trigger required \(E_{\mathrm {T}}\) above 60 GeV. A second trigger with a higher threshold of 120 GeV but looser identification requirements otherwise was used. For the 2016 data set, the main difference is that the second trigger had a threshold of 140 GeV. The trigger efficiency is at least 95% for electrons above threshold for the kinematic region used in this work.

In the muon channel, the same triggers were used for the \(\mu \mu \)jj and the \(\mu \nu \)jj channels. Events were selected by either of two single-muon triggers. For the 2015 data, one trigger required a transverse momentum of at least 26 GeV and applied an isolation criterion. To maintain a high efficiency at high \(p_{\mathrm {T}}\), a second trigger with a threshold of 50 GeV and no further requirements was used. The main difference for the 2016 data set is that slightly different isolation requirements were used for the lower-\(p_{\mathrm {T}}\) trigger. The thresholds were the same as in 2015. For the offline \(p_{\mathrm {T}}\) values considered in this analysis, the trigger efficiencies have reached their plateau values, which vary between 50 and 80% in the more central detector region and are around 90% for \(|\eta |>1.05\). The variation in trigger efficiency is due to local inefficiencies and incomplete detector coverage.

Multiple pp interactions in the same or neighbouring bunch-crossings can lead to many reconstructed vertices in the beam collision of interest. The primary vertex of the event, from which the leptons are required to originate, is defined as that with the largest sum of squared transverse momenta of its associated tracks. Selected events must contain a primary vertex with at least two associated tracks with \(p_{\mathrm {T}} > 0.4\) GeV.

3.2 Signal and background simulations

Samples of simulated events with pair-produced scalar LQs with masses between 200 and 1700 GeV, in steps of 50 GeV up to 1500 GeV and thereafter 100 GeV, were generated at next-to-leading order (NLO) in QCD [43,44,45] with the MadGraph 2.4.3 [43, 46] program using MadSpin [47] for the decay of LQs. In the simulation, leptoquarks with a charge of \(\frac{1}{3}e\) are used. Accordingly, the first generation LQ decays to either a u-quark and an electron, or a d-quark and an electron neutrino. The second generation LQ decays to either a c-quark and a muon or an s-quark and a muon neutrino. The anti-leptoquarks decay into the corresponding anti-particles. The generator output was interfaced with Pythia 8.212 [48] for the event simulation beyond the hard scattering process, i.e. the parton shower, hadronisation and underlying event, collectively referred to as UEPS. The A14 set of tuned parameters (tune) [49] was used for the UEPS modelling. The NNPDF3.0 NLO [50] parton distribution function (PDF) set was used. The coupling \(\lambda _{\mathrm {LQ} \rightarrow \ell q}\), which determines the LQ lifetime and width [13] was set to \(\sqrt{4\pi \alpha }\), where \(\alpha \) is the fine-structure constant. This value corresponds to a LQ full width of about 0.2% of its mass; LQs can thus be considered to decay promptly. Samples were generated for \(\beta \) with a value of 0.5.

Events containing W or Z bosons and associated jets [51] were simulated with the Sherpa 2.2.1 generator [52]. These event samples include off-shell production of the bosons. Matrix elements were calculated in perturbative QCD for up to two partons at NLO and four partons at leading order (LO) using the Comix [53] and OpenLoops [54] matrix element generators. Merging with the Sherpa UEPS model was performed with the ME+PS@NLO method [55]. The NNPDF3.0 PDF set [56] at next-to-next-to-leading order (NNLO) was used. Simulated processes in which a Z boson decays into leptons are hereafter termed Drell–Yan processes.

For the simulation of \(t\bar{t}\) and single-top-quark production in the Wt final state and s-channel the Powheg-Box v2 [57,58,59,60] generator was used. For the UEPS modelling for the \(t\bar{t}\) samples, Pythia 8.210 was used, with the A14 tune and the NNPDF2.3 LO [61] PDF set. The \(t\bar{t}\) samples used the NNPDF3.0 set in the matrix element calculations, the single-top samples used the CT10 PDF set. Electroweak t-channel single-top-quark events were produced with the Powheg-Box v1 generator. For the generation of single-top samples, the UEPS modelling was done using Pythia 6.428 [62] with the CTEQ6L1 [63] PDF sets and Perugia 2012 [64] tune. The value of the top mass (\(m_t\)) was set to 172.5 GeV. The EvtGen v1.2.0 [65] program was used for properties of the bottom and charm hadron decays.

Diboson processes with at least one charged lepton in the final state were simulated using the Sherpa 2.1.1 generator. All diagrams with four electroweak vertices were considered. They were calculated for up to one parton at NLO and up to three partons at LO using the Comix and OpenLoops matrix element generators and merged with the Sherpa UEPS model using the ME+PS@NLO prescription. The CT10 PDF set was used in conjunction with dedicated parton shower tuning developed by the Sherpa authors. The generator cross sections, calculated up to NLO, were used.

To model the effect of multiple proton–proton interactions in the same or neighbouring bunches (pile-up), simulated inclusive proton–proton events were overlaid on each generated signal and background event. The pile-up was simulated with Pythia 8.186 using tune A2 [66] and the MSTW2008LO PDF set [67]. Simulated events were corrected using per-event weights to describe the distribution of the average number of interactions per proton bunch-crossing as observed in data.

The detector response to events in the SM background samples was evaluated with the Geant4-based detector simulation [68, 69]. A fast simulation employing a parameterisation of calorimeter response [70] and Geant4 for the other detector components was used for the signal samples. The standard ATLAS reconstruction software was used for both simulated and pp data. Furthermore, correction factors, termed scale factors, which were derived from data were applied as event weights to the simulation of the lepton trigger, reconstruction, identification, isolation, and impact parameter selection, and of b-tagging efficiencies.

4 Object definition and event pre-selection

Electrons are reconstructed from ID tracks which are matched to energy clusters found in the electromagnetic calorimeter. The reconstruction efficiency is higher than 97% for candidates with \(p_{\mathrm {T}}\) greater than 30 GeV. An object is identified as an electron following requirements made on the quality of the associated track, shower shapes exploiting the longitudinal segmentation of the electromagnetic calorimeter, leakage into the hadronic calorimeter, the quality of the track-to-cluster matching and measurements of transition radiation made with the TRT [71]. The transverse energy of the electron candidates in the eejj (\(e\nu jj\)) channel must exceed 40 (65) GeV. Only electron candidates in the pseudorapidity region \(|\eta |<2.47\) and excluding the transition region between the barrel and endcap EM calorimeters (\(1.37<|\eta |<1.52\)) are used. The electron identification uses a multivariate technique to define different working points of selection efficiency. In the dielectron channel, the ‘medium’ identification working point with an efficiency above 90% for the candidates considered in this analysis is used. To achieve better rejection against jets misidentified as electrons, a tighter selection is used in the electron–neutrino channel. For the electron–neutrino channel, the ‘tight’ identification working point is used with an efficiency of 85% or more for candidates with \(p_{\mathrm {T}}\) above 65 GeV. The electrons are required to be isolated, using both calorimeter- and track-based criteria such that the isolation efficiency is 98%. Finally, the electron candidates have to be compatible with originating from the primary vertex.

Muon tracks are reconstructed independently in the ID and the MS [72]. Tracks are required to have a minimum number of hits in each system, and must be compatible in terms of geometrical and momentum matching. Information from both the ID and MS is used in a combined fit to refine the measurement of the momentum of each muon over \(|\eta | <2.5\) [73]. The efficiency for reconstructing muons is 98%. As for the electron channels, muon candidates are required to have \(p_{\mathrm {T}}\) > 40 GeV and \(p_{\mathrm {T}}\) > 65 GeV for the \(\mu \mu \)jj and \(\mu \nu \)jj searches respectively. They are further required to be compatible with originating from the primary vertex. A track-based isolation requirement is applied to the muons, yielding a selection efficiency of 99%.

Jets are reconstructed using the anti-\(k_{t}\) algorithm [74] with a radius parameter \(R=0.4\) from topological clusters of calorimeter cells which are noise-suppressed and calibrated to the electromagnetic scale. They are calibrated using energy- and \(\eta \)-dependent correction factors derived from simulation and with residual corrections from in situ measurements [75]. The jets must satisfy \(p_{\mathrm {T}}~>~60\) GeV and \(|\eta |~<~2.5\). Additional jet quality criteria are also applied to remove fake jets caused by detector effects [76]. Furthermore, to eliminate jets containing a large energy contribution from pile-up, jets are tested for compatibility with the hard scatter vertex with a jet vertex tagger discriminant, utilising information from the ID tracks associated with the jet [77].

Jets containing b-hadrons (b-jets) are identified with an algorithm which is based on multivariate techniques. The algorithm combines information from the impact parameters of displaced tracks and from topological properties of secondary and tertiary decay vertices reconstructed within the jet. A working point at which the b-tagging efficiency is around 77% for jets originating from a b-quark is chosen, as determined with a MC simulation of \(t\bar{t}\) processes [78, 79]. Scale factors needed to match the MC performance to data are compatible with unity [80].

The missing transverse momentum is defined as the negative vector transverse momentum sum of the reconstructed and calibrated physics objects, plus an additional soft term [81, 82]. The soft term is built from tracks that are not associated with any reconstructed electron, muon or jet, but which are associated with the primary vertex. The absolute value of the missing transverse momentum is denoted by \(E_{\mathrm {T}}^{\mathrm {miss}}\). In the \(\ell \nu jj\) channel, an \(E_{\mathrm {T}}^{\mathrm {miss}}\)-related quantity, calculated as \(S=E_{\mathrm {T}}^{\mathrm {miss}}/\sqrt{p_\mathrm {T}^{j1}+p_\mathrm {T}^{j2}+p_{\mathrm {T}}^{\ell }}\) is used to further reduce backgrounds from objects wrongly reconstructed as leptons. Here, \(p_\mathrm {T}^{j1}\) (\(p_\mathrm {T}^{j2}\)) is the \(p_{\mathrm {T}}\) of the leading (subleading) jet and \(p_{\mathrm {T}}^{\ell }\) is the \(p_{\mathrm {T}}\) of the charged lepton.

Ambiguities in the object identification which arise during reconstruction, i.e. when a reconstructed object can match multiple object hypotheses (electron, muon, jet), are resolved in several steps. First, electrons are removed if they share a track with a muon. An algorithm is also used to remove jet–lepton ambiguities based on the proximity of identified leptons and jets.

For the dilepton channel, events are selected such that they contain at least two jets and two same-flavour leptons fulfilling the requirements above for the respective flavour. Selected lepton-neutrino events have to contain at least two jets, one lepton, and a missing transverse energy of more than 40 GeV.

5 Control regions and background estimation

The major SM background processes to the LQ signal correspond to the production of (off-shell) \(Z/\gamma ^*{+}\)jets, \(W{+}\)jets and \(t\bar{t}\) events in which at least one top quark decays leptonically. These backgrounds are determined using data in selected CRs, as described below. Subdominant contributions arising from diboson, single-top, \(W\rightarrow \tau \nu \) and \(Z\rightarrow \tau \tau \) production are estimated entirely from simulation. There are also background contributions due to objects being misidentified as leptons or non-prompt leptons that are produced in the decay of hadrons. These backgrounds are collectively referred to as the fake (lepton) background and are estimated in a data-driven way where relevant.

Three CRs, independent of the SRs used for the leptoquark search (see Sect. 7), in which any high-mass LQ signal contributions are expected to be negligible are defined in order to estimate the main backgrounds. First, the \(Z/\gamma ^*{+}\)jets control region (Z CR) in the dilepton channels is defined by restricting the dilepton invariant mass, \(m_{\mathrm {\ell \ell }}\) , to values between 70 and 110 GeV. These CRs contain samples of \(Z/\gamma ^*{+}\)jets events of purity 92% (94%) for the electron (muon) channel. Next, the W control region (W CR) in the lepton–neutrino channels requires the transverse massFootnote 2, \(m_{\mathrm {T}}\) , to be between 40 and 130 GeV. In addition, \(E_{\mathrm {T}}^{\mathrm {miss}}\) has to be greater than 40 GeV and the observable S greater than 4. Events that contain one or more b-jets are vetoed to reduce the \(t\bar{t}\) contribution. The purity of the W CR is 74% (80%) for the electron (muon) channel.

To improve the description of data, a reweighting of the \(Z/\gamma ^*{+}\)jets and \(W{+}\)jets predictions by Sherpa 2.2.1 is performed. The weights in this procedure are parameterised as a function of the dijet mass, \(m_{jj}\), in the respective CR, and applied to both CRs and SRs (see Sect. 7). The parameterisation is derived by fitting the ratio between data and prediction with second-degree polynomials. Figure 2 shows distributions of \(m_{jj}\) and leading jet \(p_{\mathrm {T}}\) in the Z CR for the \(\mu \mu \)jj channel before and after reweighting. The data are better described by the Monte Carlo simulations following this reweighting. The total uncertainty (statistical, systematic and theory) is shown. The total uncertainty for CR distributions is dominated by the uncertainty in the theoretical predictions which is correlated between bins. Experimental uncertainties from jet energy scale and jet energy resolution are typically around 1%. Similarly, the Sherpa 2.2.1 predictions of \(W{+}\)jets require reweighting to match the data. The functional form of the reweighting algorithm is similar to that used for the \(Z/\gamma ^*{+}\)jets sample. The signal regions differ from the control regions mainly by the cut on \(m_{\mathrm {\ell \ell }}\) or \(m_{\mathrm {T}}\) . It is verified that applying weights derived in one mass region to the MC simulation in a different mass region provides a good description of the data in that region, indicating that the same correction is also valid in the signal region. All mass regions used in this test do not overlap with the signal regions. A systematic uncertainty is included to account for the reweighting procedure (Sect. 8).

In the statistical method used to evaluate the exclusion limits which is described in Sect. 9, the normalisations of background estimates are adjusted in the fitting procedure. These normalisations are not applied to any of the distributions shown prior to that section. Finally, the \(t\bar{t}\) control region (\(t\bar{t}\) CR) for the \(\mu \nu {} \textit{jj} \) and \(\textit{e}\nu {} \textit{jj} \) channels is defined with the same requirements on \(m_{\mathrm {T}}\) and \(E_{\mathrm {T}}^{\mathrm {miss}}\) as the W CR, but here the presence of at least two b-jets is required. The purity of the \(t\bar{t}\) CR is 86% (89%) for the electron (muon) channel. The CRs are summarised in Table 1. The Z CRs and other measurement regions are used to extract particle-level differential cross-sections (see Sect. 6).

Fig. 2
figure 2

Distributions of \(m_{jj}\) (left) and leading jet \(p_{\mathrm {T}}\) (right) in the Z CR for the \(\mu \mu \)jj channel. Data are shown together with MC contributions corresponding to \(Z \rightarrow \mu \mu \), \(t\bar{t}\), diboson (VV) and single-top (tW) processes. The MC distributions are cumulatively stacked. The bottom panels show the ratio of data to expected background. The grey hatched band represents the total uncertainty. The Sherpa 2.2.1 predictions for \(Z \rightarrow \mu \mu \) are shown before (top) and after (bottom) the reweighting procedure (see text)

Table 1 Definition of the CRs used in the dilepton and lepton–neutrino channel, respectively. Selections common to both channels are shown at the top

This analysis makes use of the matrix method [83] to estimate the background from misidentified electrons for the first-generation LQ search. This background is negligible for the second generation search. The matrix method is based on the estimation of the probabilities for prompt electrons and fake electron candidates that pass identification and isolation selections that are looser than the nominal selection described in Sect. 4, to also pass the nominal selection. These probabilities are referred to as the prompt rate and fake rate, respectively. In the loose selection, several criteria that are used to suppress fake contributions in the nominal selection are relaxed. In particular, no isolation is required in the loose selection. The prompt rate is estimated to a good approximation from MC simulations, while the fake rate is estimated from a data sample enriched in fake backgrounds. To suppress contributions from prompt electrons in this sample, events are rejected if they contain at least two electron candidates that pass the ‘medium’ identification criteria. In both channels, events are also rejected if there are two electron candidates that fulfil the ‘loose’ object definition specific to the respective channel and have an invariant mass in an interval of \(\pm 20\) GeV around the Z boson mass. The contribution from the fake electron background is at most 24%, but generally much less.

Backgrounds with non-prompt muons arise from decays of heavy-flavour hadrons inside jets. Their contribution to the dimuon final states is negligible as in the search documented in Ref. [84]. In lepton–neutrino final states, a contamination of about 5% was found [83], corresponding to about half the size of the fake electron background. Tests showed that a contamination of 5% would not alter the results of this analysis, and therefore this background is neglected. The uncertainty in the total background estimation due to neglecting this background is taken to be half the size of the background from fake electrons.

In Fig. 3, the distributions of \(m_{\mathrm {LQ}}^{\mathrm {min}}\) and the dilepton invariant mass in the Z CR are shown for the \(\mu \mu \)jj and eejj channels. The quantities \(m_{\mathrm {LQ}}^{\mathrm {min}}\) ane \(m_{\mathrm {LQ}}^{\mathrm {max}}\) are defined as the lower and higher of the two invariant masses which can be reconstructed using the two lepton–jet pairs in the dilepton channels. In the analysis, the lepton–jet pairing is chosen such that the absolute mass difference between the two \(\mathrm {LQ}\) candidates is minimised. The full set of experimental and theoretical uncertainties is considered for the uncertainty band. The spectra correspond to the predictions prior to the fit described in Sect. 9.

Distributions for the e \(\nu \)jj and \(\mu \nu \)jj channels are similarly shown for the W and \(t\bar{t}\) CRs in Figs. 4 and  5, respectively. The \(E_{\mathrm {T}}^{\mathrm {miss}}\) and \(m_{\mathrm {T}}\) spectra are shown in Fig. 4, while the \(m_{jj}\) and \(m_{\mathrm {T}}\) spectra are shown in Fig. 5. The data and the predictions in the various CRs are in agreement within uncertainties (described in Sect. 8) following the reweighting of the Sherpa 2.2.1 predictions for \(Z/\gamma ^*{+}\)jets and \(W{+}\)jets processes. After the reweighting, the overall normalisation of the simulation is in very good agreement with the data in the \(V{+}\)jets CRs.

Fig. 3
figure 3

Distributions of \(m_{\mathrm {LQ}}^{\mathrm {min}}\) and \(m_{\mathrm {\ell \ell }}\) in the Z CR for the eejj (top) and \(\mu \mu \)jj (bottom) channels. Data are shown together with background contributions corresponding to \(t\bar{t}\), diboson (VV), \(Z \rightarrow ee\), \(Z \rightarrow \mu \mu \), single-top (tW) processes, and fake electrons. The background distributions are cumulatively stacked. The bottom panels show the ratio of data to expected background. The grey hatched band represents the total uncertainty. A reweighting as a function of \(m_{jj}\) is applied to the \(Z{+}\)jets simulation. The other MC predictions are not reweighted

Fig. 4
figure 4

Distributions of \(E_{\mathrm {T}}^{\mathrm {miss}}\) and \(m_{\mathrm {T}}\) in the W CR for the e \(\nu \)jj (top) and \(\mu \nu \)jj (bottom) channels. Data are shown together with background contributions corresponding to \(t\bar{t}\), diboson (VV), \(Z\rightarrow ee\), \(Z\rightarrow \mu \mu \), \(W\rightarrow \mu \nu \), \(W\rightarrow e\nu \), \(W\rightarrow \tau \nu \), single-top production (tW, tq) processes, and fake electrons. The background distributions are cumulatively stacked. The bottom panels show the ratio of data to expected background. The grey hatched band represents the total uncertainty. A reweighting as a function of \(m_{jj}\) is applied to the \(W{+}\)jets and \(Z{+}\)jets simulation. The other MC predictions are not reweighted

Fig. 5
figure 5

Distributions of \(m_{jj}\) and \(m_{\mathrm {T}}\) in the \(t\bar{t}\) CR for the e \(\nu \)jj (top) and \(\mu \nu \)jj (bottom) channels. Data are shown together with background contributions corresponding to \(t\bar{t}\)  diboson (VV), \(Z\rightarrow ee\), \(Z\rightarrow \mu \mu \), \(W\rightarrow \mu \nu \), \(W\rightarrow e\nu \) \(W\rightarrow \tau \nu \), single-top production (tW, tq) processes, and fake electrons. The background distributions are cumulatively stacked. The bottom panels show the ratio of data to expected background. The grey hatched band represents the total uncertainty. A reweighting as a function of \(m_{jj}\) is applied to the \(Z{+}\)jets and \(W{+}\)jets simulations. The other MC predictions are not reweighted

6 Extraction of particle-level cross-sections from the measurement regions

Measurements of the particle-level fiducial and differential cross-sections are made in several regions related to the search CRs described in the previous section. These measurement regions (MRs) correspond to kinematic regions where new physics signals have already been excluded. They are chosen because several searches with related final states use similar selections for their CRs (see for example Ref. [38]).

6.1 Particle-level objects

In order to extract fiducial and differential cross-sections in the MRs, particle-level object selections are defined, and chosen to be as close as possible to their detector-level counterparts described in Sect. 4.

Prompt dressedFootnote 3 particle-level electrons and muons are used to define the fiducial region. They are required to have \(p_{\mathrm {T}}>40{\mathrm {\ GeV}}\) and \(|\eta |<2.5\).

Particle-level jets are formed by clustering stableFootnote 4 particles using the anti-\(k_t\) (\(R=0.4\)) algorithm [74], excluding muons and neutrinos. Particle-level jets are required to have \(p_{\mathrm {T}}>60{\mathrm {\ GeV}}\) and \(|\eta |<2.5\). The electrons used in the jet clustering are not dressed. To match the requirements of the detector-level jet-lepton ambiguity resolution procedure, any particle-level leptons within \(\Delta R <0.4\) from a selected jet are removed.

Table 2 Requirements at particle-level for all measurement regions (MRs) where particle-level differential cross-sections are extracted. The purity refers to the proportion of the particle-level yield of the MR which is accounted for by the dominant process. The variable \(S_{T}\) refers to the scalar sum of the \(p_{\mathrm {T}}\) of the jets and leptons

6.2 Measurement region definitions

The first two MRs are identical to the eejj and \(\mu \mu \)jj Z CRs described in Table 1. The third is an \(e\mu \)jj region where the leptons are required to have different flavours, and no requirement on the dilepton mass is applied. Finally, three additional regions are defined as for the eejj, \(\mu \mu \)jj and \(e\mu \)jj MRs, and labelled as extreme with the additional requirement that the scalar sum of the \(p_{\mathrm {T}}\) of the leptons and the two leading jets (\(S_{\mathrm {T}}\)) be above 600 GeV. The requirements for events entering all MRs are summarised in Table 2.

The MRs are defined both at detector level and particle level, where the requirements are identical but use the objects passing the detector-level and particle-level selections respectively. The \(m_{jj}\) reweighting is also applied to events from the simulated Drell–Yan sample at particle level, where the particle-level value of \(m_{jj}\) is used to calculate the weight. In each MR, the cross-section measurements are made for dominant processes. The contributions from all other processes are subtracted from the data yield, as explained in Sect. 6.3.

6.3 Particle-level measurement method

Each MR is designed such that it is dominated by a single SM process (accounting for more than 85% of the yield, see Table 2), for which the cross-section can be measured exclusively. A simple procedure to extract particle-level cross-sections, based on bin-by-bin correction factors, is used. In this method, the particle-level differential cross-section for the dominant process p in bin i of the distribution of variable X can be expressed as:

$$\begin{aligned} \frac{{\mathrm {d}} \sigma _i^p}{{\mathrm {d}}X} = \frac{\left( N_i - \sum _{q\ne p} R_i^q\right) \cdot \frac{T_i^p}{R_i^p}}{w_i \cdot L} , \end{aligned}$$
(1)

where the sum runs over all sub-dominant SM processes which contribute to the yield, \(N_i\) is the number of events observed in the bin, \(R_i^p\) (\(T_i^p\)) is the yield at detector level (particle level) for process p, \(w_i\) is the width of the bin and L is the integrated luminosity of the data sample. The bin-by-bin correction factors \(T_i^p/R_i^p\) may be thought of as the inverse of the reconstruction efficiency in bin i, assuming no bin-by-bin migrations. The bin-by-bin correction procedure doesn’t account for migration between bins due to resolution effects when going from particle-level to detector-level distributions. For this reason, the binning of the distributions is chosen to be much broader than the experimental resolution to minimise the migration between bins. This binning is optimised for the measurements to ensure that at least 90% of events passing the particle-level selection remain in the same bin at detector level (for events passing both selections). This binning is not used for the search.

Differential cross-section measurements are made for the following observables: \(p_{\mathrm {T}}\) of the dilepton system; \(\Delta \phi \) between the leptons; minimum \(\Delta \phi \) between lepton and leading jet; minimum \(\Delta \phi \) between lepton and subleading jet; \(S_{\mathrm {T}}\); leading jet \(p_{\mathrm {T}}\); subleading jet \(p_{\mathrm {T}}\); \(\Delta \phi \) and \(\Delta \eta \) between the leading and subleading jets; scalar sum of \(p_{\mathrm {T}}\) of leading and subleading jets (\(H_\mathrm {T}\)); and invariant mass of the dijet system.

7 Multivariate analysis and signal regions

In order to discriminate between signal and background in the LQ search, the TMVA [35] implementation of a BDT is used. A number of variables are expected to provide discrimination between signal and background. Generally, the \(p_{\mathrm {T}}\) of the leptons and jets as well as related variables such as \(S_{\mathrm {T}}\) will possess higher values for the signal than for the background, in particular for high LQ mass. Since the signature arises from the decay of parent LQs with well-defined masses, mass-sensitive discriminating variables are further candidates for BDT input variables. In the dilepton channels, the lepton-jet masses \(m_{\mathrm {LQ}}^{\mathrm {min}}\) and \(m_{\mathrm {LQ}}^{\mathrm {max}}\) are reconstructed, as described in Sect. 5. In the lepton–neutrino channels, the mass of one \(\mathrm {LQ}\) (\(m_{\mathrm {LQ}}\)) and the transverse mass of the second \(\mathrm {LQ}\) (\(m_{\mathrm {LQ}}^{\mathrm {T}}\)) can be reconstructed. The transverse mass is defined as \(m_{\mathrm {LQ}}^{\mathrm {T}}=\sqrt{2\cdot p_T^j \cdot E_{\mathrm {T}}^{\mathrm {miss}}\cdot (1-\cos (\Delta \phi (j,E_{\mathrm {T}}^{\mathrm {miss}})))}\), where \(\Delta \phi (j,E_{\mathrm {T}}^{\mathrm {miss}})\) is the azimuthal angle between the jet and the direction of missing transverse momentum vector. The procedure is analogous to that used in the dilepton channel; both possible pairings of jet and lepton (or \(E_{\mathrm {T}}^{\mathrm {miss}}\)) are tested and the pairing that results in the smaller absolute difference between \(m_{\mathrm {LQ}}\) and \(m_{\mathrm {LQ}}^{\mathrm {T}}\) is chosen.

All background processes are used in the BDT training. For the electron and muon channels separately, the training is done for each \(\mathrm {LQ}\) mass hypothesis for which signal simulations were produced. To ensure independence from the CRs, in the dilepton channel, only events with a dilepton invariant mass above 130 GeV are considered in the training. Similarly, only events with a transverse mass greater than 130 GeV are considered in the lepton–neutrino channel. Here, to further suppress the fake-electron contribution, it is required that \(E_{\mathrm {T}}^{\mathrm {miss}}\)\(>150\) GeV and that \(S>3\). These requirements define the training regions (TRs), i.e. event samples used to train the BDTs.

To ensure that no bias is introduced, two BDTs are used for each mass point; one is trained on one half of the events and evaluated on the other half, and vice versa for the second BDT. Different methods of dividing the simulation samples into training and test samples are tested and compared, giving consistent results such that no overtraining or other biases are introduced in the process.

Table 3 lists the variables that were found to give optimal sensitivity, ordered by their discriminating power in the BDT. In the dilepton channel, these are both reconstructed \(\mathrm {LQ}\) candidate masses, \(m_{\mathrm {LQ}}^{\mathrm {min}}\) and \(m_{\mathrm {LQ}}^{\mathrm {max}}\), the dilepton invariant mass, \(m_{\mathrm {\ell \ell }}\) , the subleading jet \(p_{\mathrm {T}}\) (\(p_\mathrm {T}^{j2}\)) and the subleading lepton \(p_{\mathrm {T}}\) (\(p_\mathrm {T}^{\ell 2}\)). The lepton–neutrino channel also makes use of both of the LQ mass variables, \(m_{\mathrm {LQ}}\) and \(m_{\mathrm {LQ}}^{\mathrm {T}}\), and the other variables are the transverse mass, \(m_{\mathrm {T}}\) , the missing transverse momentum, \(E_{\mathrm {T}}^{\mathrm {miss}}\), the subleading jet \(p_{\mathrm {T}}\), and the charged-lepton \(p_{\mathrm {T}}\), \(p_{\mathrm {T}}^{\ell }\). Distributions of a selection of the variables used for training are shown in Appendix B. It was checked and confirmed that not only the individual input variables but also their pairwise correlations are modelled well in the simulation.

Table 3 BDT input variables in the dilepton and lepton–neutrino channels

Output distributions for the BDT discriminant in the TRs are shown in Fig. 6 for a signal sample of LQ mass 1.3 TeV. The background distributions are shown before the fit (see Sect. 9). The contribution of a LQ of mass 1.3 TeV, normalised to the production cross-section, is also shown for the four channels. The background includes the reweighted predictions of \(Z/\gamma ^*{+}\)jets and \(W{+}\)jets. An increased separation between signal and background at higher values of the BDT output is seen.

Fig. 6
figure 6

Output distributions for the BDT in the training regions (TRs) for the four search channels: a eejj, b e \(\nu \)jj, c \(\mu \mu \)jj and d \(\mu \nu \)jj. Data are shown together with pre-fit background contributions corresponding to \(t\bar{t}\), diboson (VV), \(Z\rightarrow ee\), \(Z\rightarrow \mu \mu \), \(W\rightarrow \mu \nu \), \(W\rightarrow e\nu \), \(W\rightarrow \tau \nu \), single-top production (tW, tq) processes, and fake electrons. The prediction for a LQ with mass 1.3 TeV is also shown. The predicted background distributions are cumulatively stacked. The bottom panels show the ratio of data to expected background. The grey hatched band represents the total uncertainty

The SR phase-space is a subset of that considered in the TRs. The same object and event selections as in the TRs are made in the SRs, but in addition a requirement is placed on the BDT score. As stated above, each BDT is trained on half of the events in the TR and evaluated on the other half, and the obtained BDT score distribution is used to determine the SR. In this way, the SR determined based on a given BDT does not overlap with the TR used for that BDT. For the dilepton and lepton–neutrino channels separately, the SRs are defined as a single bin in the output score distribution of the respective BDT. The bin range is determined for each mass point separately by maximising the quantity [85] \({Z=\sqrt{(s+b)\log (1+s/b)-s}}\), where s and b are the signal and background expectations, respectively. In this optimisation, a number of further requirements are placed on the background expectations: there have to be at least two background events in the chosen region; the statistical uncertainty of the number of background events estimated by the simulation has to be less than 20%; if less than ten background events are expected, the statistical uncertainty is required to be less than 10%. This is done to avoid selecting bins with very low background expectation and thus to increase the fit stability.

In total, there are 58 SRs in the electron and the muon channel each: two SRs (\(\ell \ell jj\) and \(\ell \nu jj\)) for each of the 29 mass hypotheses considered. The fraction of signal events that is left after acceptance selections and losses due to detector inefficiencies, is estimated from the simulation to rise from less than 1% at low masses (around 200 GeV) to 60–70% for masses around 1500 GeV.

8 Systematic uncertainties

The systematic uncertainties for the search are divided into two categories: experimental uncertainties and theoretical uncertainties in the background and the signal. Most of the sources of uncertainties are common to both the search and the measurement. In addition to the common set of uncertainties, the measurement incorporates MC modelling uncertainties from closure tests. The search applies uncertainties to normalisation factors of the major backgrounds and their extrapolation to the SRs; these are evaluated in the CR-only profile likelihood fit in which the normalisation factors are determined. In Sects. 8.1 and 8.2, the sources of experimental and theoretical uncertainies are outlined, respectively. The impact of the major uncertainties on the search and the measurements are described in Sects. 8.3 and 8.4, respectively.

8.1 Experimental uncertainties

The uncertainty in the combined 2015 + 2016 integrated luminosity is 2.1%. It is derived, following a methodology similar to that detailed in Ref. [86], and using the LUCID-2 detector for the baseline luminosity measurements [87], from calibration of the luminosity scale using x–y beam-separation scans.

An uncertainty due to the reweighting of the vertex multiplicity to describe pile-up is estimated by scaling the average number of interactions per bunch-crossing in data by different factors.

The major sources of uncertainties affecting the electron measurement derive from the electron energy scale and resolution measurements. The analysis also considers uncertainties due to the modelling of the efficiencies of the four electron selection criteria: trigger, reconstruction, identification and isolation.

Sources of experimental uncertainties related to muon reconstruction and calibration include uncertainties in the determination of the MS momentum scale, MS momentum resolution, ID momentum resolution and additional charge-dependent corrections. Furthermore, uncertainties in the determination of the four efficiency scale factors, introduced at the end of Sect. 3.2, for trigger, identification, isolation and track-to-vertex matching are taken into account.

There are two main sources of uncertainty related to the jet reconstruction: jet energy scale (JES) and jet energy resolution (JER). Another source of jet-related uncertainty corresponds to corrections made for the b-tagging efficiency. The uncertainty due to the JES is largely derived from various in situ techniques. Four components of this uncertainty are considered.

The uncertainties due to the electron, muon and jet energy scale and resolution are propagated to the estimation of \(E_{\mathrm {T}}^{\mathrm {miss}}\). In addition, there is an uncertainty from the soft term, which has three components: one corresponding to the \(E_{\mathrm {T}}^{\mathrm {miss}}\) scale and two components for the \(E_{\mathrm {T}}^{\mathrm {miss}}\) resolution uncertainty. The resolution uncertainty is split into components parallel and perpendicular to an axis in the transverse plane which is defined along the direction of a vectorial sum of all the hard objects in the event (electrons, muons and jets).

An uncertainty in the fake rate of electrons is determined by varying the MC prediction that is used to remove the contributions of prompt electrons from the fake-enriched sample by 30%. The effect is propagated to the fake background estimate. The uncertainty due to neglecting the muon fake background is taken to be half the size of the electron fake background, see Sect. 5.

8.2 Theoretical uncertainties

The Sherpa 2.2.1 \(V{+}\)jets samples include event weights reflecting variations of the nominal PDF set (NNPDF3.0) and the use of two different PDF sets: MMHT2014NNLO68CL [88] and CT14NNLO [89]. The NNPDF intra-PDF uncertainty is estimated as the standard deviation of the set of 101 NNPDF3.0 eigenvectors. This procedure takes into account the effect of varying \(\alpha _{\text {S}} \) by \(\pm \,0.001\) around its nominal value of 0.119. The envelope of the differences between the nominal NNPDF set and the other two PDF sets is used as an additional uncertainty.

The samples include weights for seven sets of variation of the renormalisation (\(\mu _\textsc {r} \)) and factorisation (\(\mu _\textsc {f} \)) scale, i.e. varying the scales up and down by a factor of two, either together or independently. The envelope of all these variations is taken as an estimate of the scale uncertainty.

A further uncertainty from the reweighting of the \(V{+}\)jets simulations in \(m_{jj}\) is considered. The full difference between the reweighted and the unweighted distributions is taken as an estimate of the uncertainty induced by the reweighting.

The uncertainties in the modelling of the production of \(t\bar{t}\) are assessed from a number of alternative simulation samples. Differences between different generators and different models for fragmentation and hadronisation are considered. In addition, parameters affecting initial- and final-state radiation are also varied. For the \(t\bar{t}\) sample the default PDF set is NNPDF3.0 and the PDF uncertainty is estimated from the PDF4LHC15 prescription [90].

8.3 Uncertainties for the search

As described in Sect. 9 three CRs are fit simultaneously to obtain normalisation factors. Sources of theoretical uncertainty described above are taken into account in this procedure. Furthermore, theory uncertainties in the SRs correspond to variations in the output BDT spectra and are thus taken into account in the final limit setting calculation.

The electron energy scale and reconstruction affect the total background yield in the SRs for the electron channels by 2–14%. Similarly uncertainties due to the muon reconstruction and calibration lead to background uncertainties of between 2 and 20% for the muon channels. The impact of uncertainties related to the JES (JER) lead to uncertainties on the predicted background of 5–20% (2–20%). The effect of the soft-term uncertainties in the \(E_{\mathrm {T}}^{\mathrm {miss}}\) estimations leads to variations of the background estimate from less than 1% to 15%. Generally, the experimental uncertainties in the signal yields in the SRs are not larger than 2%.

Uncertainties of approximately 100% for the fake electron background are estimated at low values of lepton-jet mass in the TR and about 20% at masses above 1 TeV. The statistical uncertainty in the electron fake rate determination is below 10%. As discussed earlier, the contribution from fake muons is typically half the size of the electron fake background, and this is used as a (conservative) estimate of the uncertainty in the total background due to neglecting the muon fake contribution.

The final modelling uncertainties in the background yields for \(Z{+}\)jets and \(W{+}\)jets are around 2–10%, with the scale uncertainty typically being the dominant source. Signal predictions have an uncertainty of a similar size. The background from \(t\bar{t}\) production is affected by a theoretical unceratinty of up to \(50\%\), the large value arising due to the statistical imprecision of the MC samples.

8.4 Uncertainties for the cross-section measurements

The uncertainties described below are those which apply to the particle-level cross-section measurements in the MRs, in addition to the relevant uncertainties from the previous sections. As is described in Sect. 6, the cross-section calculation involves the subtraction of non-dominant processes from the observed data, and the application of bin-by-bin correction factors. The experimental uncertainties play a role in both of these steps, since they affect the detector-level yield of both the non-dominant processes (when subtracting from data) and the dominant process (when calculating the correction factors). The dominant experimental uncertainties are the JES and the lepton efficiencies, which have effects of up to 5% and 2% on the fiducial cross-sections respectively. Theory and modelling uncertainties in the dominant process affect both the particle-level and detector-level yields, and therefore largely cancel out when calculating the correction factors T / R (see Eq. 1).

The theory uncertainties in the dominant processes in the \(Z{+}\)jets and \(t\bar{t}\) MRs therefore have only a small impact (\(\le 2\%\)) on the fiducial cross-sections. The theory uncertainties in the estimation of non-dominant processes generally have a minimal impact since these processes usually account for only small fractions of the MR yield. The exception is the case of the \(t\bar{t}\) background to the Z MRs, which also contribute a small but non-negligible theory uncertainty.

In addition, bin-by-bin correction uncertainties are derived by comparing the nominal correction factors with those obtained using re-weighted \(V{+}\)jets and alternative \(t\bar{t}\) samples. The size of these modelling uncertainties ranges between 0.5 and \(4\%\) of the fiducial cross-sections.

The overall uncertainty of the fiducial cross-sections, accounting for all sources of systematic uncertainty, is between 5% and 9% depending on the MR.

9 Statistical procedure for the search

The results of the analysis are interpreted using a profile likelihood method as implemented in HistFitter [91], in particular using the asymptotic formulae from Ref. [85]. The signal and backgrounds are described by a binned probability density function (p.d.f.) built using either the three CRs described in Sect. 5, or the two SRs, one in the dilepton channel and one in the lepton–neutrino channel, as described in Sect. 7. All CRs and SRs consist of only a single bin.

The CRs are enriched in a specific background process and have a negligible signal contamination; they are thus used to normalise the predicted backgrounds to data, whereas the SRs drive the signal extraction. The p.d.f. for the fit to the SRs includes the signal strength, i.e. the scaling factor with respect to the predictions for the signal cross-section, as the parameter of interest. Systematic uncertainties are incorporated in the p.d.f. as statistical nuisance parameters. They are introduced as shape uncertainties, i.e. they can only affect the relative size of a given background in different regions while preserving the overall number of expected events. This is done in order to remove strong correlations with the normalisation factors for the main backgrounds that are extracted from the CRs and which are used to scale MC yields to match data. Experimental uncertainties are treated as being fully correlated between different physics processes and fit regions. Theoretical or modelling uncertainties are treated as correlated in different fit regions, but uncorrelated between processes. The statistical uncertainties of the MC samples are included in the p.d.f. as nuisance parameters by constraint terms described by a \(\Gamma \) distribution. The fit is performed in two stages: first, only the CRs are included in the fit in order to extract the normalisation factors for the three main backgrounds. These normalisation factors are then applied to the main backgrounds in the second step, which is a fit to only the SRs. The uncertainty in the normalisation factors from the CR fit is introduced in the form of Gaussian nuisance parameters. The best-fit central values and constraints on other nuisance parameters are not transferred from the CRs to the SRs. The optimal value and error of the signal strength and nuisance parameters as well as their correlations are determined simultaneously when the p.d.f. is fitted to data. This conservative two-step approach is chosen for its simplicity compared with the case of a simultaneous fit to all regions.

For some systematic uncertainties, certain adjustments are made before they are incorporated in the fit, as described in the following. Uncertainties larger than 50% are capped at 50%, because their large size is due to the statistical uncertainty of simulation samples used to estimate them. This applies in particular to the theoretical uncertainties in the \(t\bar{t}\) estimation. Uncertainties with very asymmetric up and down variations are symmetrised, using the larger of the two variations. In the SR fit, experimental uncertainties that before normalisation have an effect smaller than 2% on the total background or the signal expectation are not considered in the fit for the background or the signal, respectively. Other nuisance parameters are not considered for certain background processes if their effect is less than 2% on that process before normalising the uncertainty.

10 Results

10.1 Results of the cross-section measurements

The measured fiducial particle-level cross-sections are shown in Table 4 for each of the six MRs. The values in all MRs are found to agree with the generator predictions within uncertainties. The uncertainties in the measured cross-sections are in all cases dominated by experimental uncertainties (in particular the jet energy scale). The theory uncertainties in the generator predictions are dominated by variations in the cross-sections, arising from QCD scale variations in the phase-space regions of the measurements. In both the standard and extreme MRs, the measured fiducial cross-sections in the Z MRs are found to agree between electron and muon decay modes, as expected.

Examples of the particle-level differential cross-sections are shown in Figs. 7 and 8, where they are compared with the cross-sections predicted from the nominal simulation without the \(m_{jj}\) reweighting. Predictions from alternative generators are also shown in each case. Differential cross-sections for a larger list of variables are included in Appendix A and in the HEPData [92] record for this paper.

In the eejj and \(\mu \mu \)jj MRs (and their extreme counterparts) the predicted differential cross-sections for variables involving jet energies exhibit a degree of mis-modelling by the Sherpa 2.2.1 simulation. This can be seen in Figs. 7 and 8, for the leading jet \(p_{\mathrm {T}}\) as an example. Similar mis-modelling is seen for \(m_{jj}\), \(H_\mathrm {T}\), subleading jet \(p_{\mathrm {T}}\), and \(S_\mathrm {T}\) differential measurements. All other quantities are shown to be well modelled by the nominal Sherpa 2.2.1 simulation. In all cases, the prediction of an alternative generator, MadGraph5_aMC@NLO (with events generated at NLO accuracy for \(Z+0, 1, 2 \) jets, showered using Pythia8 with the A14 tune, and for which different jet multiplicities were merged using the FxFx prescription [95]), is also shown. Typically, the MadGraph5_aMC@NLO prediction is also discrepant with the data for the same observables, but often in a direction opposite to that of the Sherpa 2.2.1 prediction. This suggests that the discrepancies are caused by choices of the generator parameters, which would have been tuned in a different phase-space region given by previous \(Z+\text {jets}\) measurements. The measurements provided by this paper will therefore help to validate new tunes of the generators to ensure better agreement in this region for future searches. These disagreements in modelling are covered by the reweighting uncertainty when evaluating limits for the search.

For the \(e\mu \)jj MRs, the predictions from the nominal Powheg simulation are found to agree with the data for all observables aside from the dilepton \(p_{\mathrm {T}}\), where a slight mis-modelling is observed in the high tail, as can be seen in the examples shown in Figs. 7 and 8. The predictions of an alternative generator (MadGraph5_aMC@NLO) are also shown, and these do agree with the data for this observable. This points to the generator tuning choice as the source of the disagreement.

Table 4 Measured fiducial cross-section for the dominant process in each of the MRs, with uncertainties broken down into the statistical, experimental, and theoretical components. The generator predictions with their uncertainties are obtained from Sherpa 2.2.1 for the \(Z\rightarrow ee\) and \(Z\rightarrow \mu \mu \) samples, and Powheg+Pythia8 for the \(t\bar{t}\) sample. Alternative generator predictions are also provided, obtained using MadGraph5 FxFx and MadGraph5+Pythia8 for the Z MRs and \(t\bar{t}\) MRs, respectively
Fig. 7
figure 7

Examples of the measured particle-level differential cross-sections in the eejj, \(\mu \mu \)jj and \(e\mu \)jj measurement regions, exclusively for the dominant process in each channel (\(Z\rightarrow ee\), \(Z\rightarrow \mu \mu \) and \(t\bar{t}\), respectively). The MC prediction for the dominant process is also shown, with no \(m_{jj}\) reweighting applied. The red band represents the statistical component of the total uncertainty on the measurement which is indicated by the error bar on each point. The variable \( min \Delta \Phi (j_{0},l) \) refers to the minimal difference in \(\phi \) between the leading jet and a prompt lepton

Fig. 8
figure 8

Examples of the measured particle-level differential cross-sections in the extreme eejj, \(\mu \mu \)jj and \(e\mu \)jj measurement regions, exclusively for the dominant process in each channel (\(Z\rightarrow ee\), \(Z\rightarrow \mu \mu \) and \(t\bar{t}\), respectively). The MC prediction for the dominant process is also shown, with no \(m_{jj}\) reweighting applied. The red band represents the statistical component of the total uncertainty on the measurement which is indicated by the error bar on each point. For the leading jet \(p_{\mathrm {T}}\) measurement, the uncertainty in the second bin is smaller because the bin-by-bin correction factor from the alternative \(t\bar{t}\) samples agrees more closely with the nominal ones than in neighbouring bins. The variable \( min \Delta \Phi (j_{0},l) \) refers to the minimal difference in \(\phi \) between the leading jet and a prompt lepton

10.2 Results of the leptoquark search

The data event yields (black points) in a representative subset of the SRs are shown in Fig. 9 together with the predicted background yields as evaluated with the fit. The bottom panel shows the significance [85] of the deviations, assuming Poisson statistics and taking only statistical uncertainties into account. Even without including systematic uncertainties, no significant excess above the SM expectation is observed in any of the SRs. The modified frequentist \({\mathrm {CL}}_{\text {s}}\) method [94] is used to set limits on the strength of the LQ signal.

Fig. 9
figure 9

Data and background yields in the SR for different values of LQ mass for a first- and b second-generation LQs. Each bin in these plots corresponds to one SR, the bin label indicating the channel, i.e. dilepton or lepton–neutrino. Two consecutive bins show the two SRs for a given mass hypothesis, that is indicated at the bottom of the plots. Mass points from 375 GeV to 1.5 TeV are shown. The background contributions correspond to  \(t\bar{t}\), \(W{+}\)jets, diboson (VV), \(Z{+}\)jets, single-top (WtWq) processes, and fake electrons. The background distributions are cumulatively stacked. The grey band indicates the total uncertainty in the background estimate after the fit. The bottom panel shows the significance of the deviations, taking only statistical uncertainties into account

Each BDT uses as inputs variables that reconstruct the LQ mass. Since these have a high discriminating power, the sensitivity of the SRs exhibits a strong dependence on the LQ mass below values of \(m_{\mathrm {LQ}}\) around 600 GeV, where the intervals between simulated mass points are greater than the typical mass resolution for adjacent masses (the converse is true for masses above around 600 GeV). Therefore, in this regime a simple interpolation of acceptance between simulated mass points does not give a reliable estimate of the limit in the mass interval. To account for this limitation, results are presented separately for low and high mass regions, defined as being less than and greater than 600 GeV, respectively.

A coarse scan of the mass is performed in the low mass region in intervals of 50 GeV. Exclusion limits that are obtained at these scan points are used as conservative limits for the intermediate masses between these points. While for masses above \(\sim 600\) GeV the mass-scan intervals are smaller than the mass resolution at the test points, the intervals used below \(\sim 600\) GeV are larger than the mass resolution. Therefore, above \(\sim 600\) GeV a simple interpolation is used to obtain the limits in the mass intervals. In the range below \(\sim 600\) GeV this is not appropriate and a different approach is used. Since the acceptance of a given SR is highest for the mass point which the SR is designed for and lower at the neighbouring mass points, the limit is estimated for a given mass interval, from \(m_1\) to \(m_2\), in the following way; two limits are evaluated for \(m_1\) and \(m_2\) using the SR defined by \(m_1\) (i.e. the respective BDT trained for the mass hypothesis \(m_1\)). Similarly, two estimations are calculated for the SR defined by \(m_2\). For each SR, the weaker of the two limits is retained, leaving two limits to consider, one from each SR. Of these two remaining results from the two SRs, the stronger limit is used.

Table 5 Expected and observed 95% CL lower limits on first- and second-generation LQ masses in units of GeV for different values of the branching ratio \(\beta \) into a charged lepton and quark
Fig. 10
figure 10

Upper limits (observed and expected) on the pair-production cross-section for (top) first- and (bottom) second-generation LQs normalised to the predicted cross-section (\(\sigma /\sigma _\text {th}\)), as a function of LQ mass for an assumed value of \(\beta =0.5\). Limits are presented for low (left) and high (right) mass regions, which correspond to masses less than and greater than 600 GeV, respectively. The \(\pm ~1\sigma \) (\(\pm 2\sigma \)) uncertainty bands on the expected limit represent all sources of systematic and statistical uncertainty

Fig. 11
figure 11

Upper limits (observed and expected) on the branching ratio for (top) first- and (bottom) second-generation LQs into a lepton and quark as a function of LQ mass. Limits are presented for low (left) and high (right) mass regions, which correspond to masses less than and greater than 600 GeV, respectively. The \(\pm ~1\sigma \) (\(\pm 2\sigma \)) uncertainty bands on the expected limit represent all sources of systematic and statistical uncertainty

The 95% CL upper limits on the cross-section for scalar-LQ pair production, normalised to the predicted cross-section, are presented as a function of \(m_{\mathrm {LQ}}\) in Fig. 10 for first- and second-generation LQs for an assumed value of \(\beta =0.5\). Expected limits with their one and two standard-deviation bands are also shown. The observed limits are consistent with the expected limits for all mass points for both channels. Lower mass limits of about 1.25 TeV are obtained, representing an increase of around 400 GeV compared with earlier work [21].

Exclusion contours in the \(\beta \)\(m_{\mathrm {LQ}}\) plane are shown in Fig. 11 for different values of \(\beta \) for the first two LQ generations. The sensitivity in \(\beta \) is greatest for the lowest mass region considered (\(\sim \) 200 GeV) for which the search is sensitive to \(\beta \) values around \(10^{-2}\). For \(\beta =1\), a mass of 1400 GeV (1560 GeV) can be excluded for first-generation (second-generation) LQs. Limit values are also given for first- and second-generation LQs in Table 5. The considerably stronger observed limit for the second generation is due to a downward fluctuation in the data in the relevant SRs.

11 Summary and conclusions

Searches for pair production of first- and second-generation scalar leptoquarks in pp collisions at \(\sqrt{s} = 13\) TeV have been made using the ATLAS detector at the LHC. The searches exploit a data set corresponding to 36.1 fb\(^{-1}\) of integrated luminosity and probe the lepton–quark and lepton–neutrino LQ decay channels. No significant excess above the SM background expectation is observed in any channel and exclusion limits have been evaluated. The results presented here significantly extend the sensitivity in mass compared to previous ATLAS results, and yield an exclusion very similar to that found by the CMS experiment using a dataset of a similar size [27, 28]. Within the minimal Buchmüller–Rückl–Wyler model and assuming a branching ratio for the decay into a charged lepton and a quark of 50%, leptoquarks with masses up to 1.29 TeV are excluded at 95% CL for first generation leptoquarks, and up to 1.23 TeV for second generation leptoquarks.

In addition to the search, measurements have been made of particle-level fiducial and differential cross-sections for six measurement regions related to the control regions used for the search. Two measurement regions are identical to the Z control regions used in the search, and their fiducial cross-sections are measured to be \(3.28\pm 0.22\) pb and \(3.32\pm 0.23\) pb for the \(Z\rightarrow ee\) and \(Z\rightarrow \mu \mu \) channels respectively. The measurements agree with the cross-sections predicted by Sherpa 2.2.1, which are \(3.24 \pm 1.02\) pb and \(3.12 \pm 1.01\) pb respectively. A third measurement region is dominated by \(t\bar{t}\rightarrow e\mu \), where the fiducial cross-section is measured to be \(1.50\pm 0.13\) pb, compared to \(1.60 \pm 0.38\) pb predicted by Powheg+Pythia8. Measurements are also made in three other regions where the scalar sum of the \(p_{\mathrm {T}}\) of the jets and leptons is above \(600{\mathrm {\ GeV}}\). These so-called “extreme” regions may be useful for generator tuning since they represent regions which are utilised in searches but where measurements are rarely made.

In addition to the inclusive cross-section measurements, differential measurements are made for eleven variables: the transverse momenta of the leading and subleading jets, the minimum angles between the leading and subleading jets with a lepton, the dilepton transverse momentum, the opening angle between the leptons, the scalar sum of the jet transverse momenta, the scalar sum of the jet and lepton transverse momenta, and the opening angles in \(\eta \) and \(\phi \) between the two leading jets. In the Z measurement regions, the differential cross-sections for variables involving leptons are typically found to be well-modelled by the nominal generators and the alternative generators shown typically give reasonable but weaker agreement. The measurements involving jet energies and angles are found to exhibit a degree of mis-modelling by both the nominal and alternative simulations. The nominal and alternative generator predictions are often discrepant with respect to the data in different directions, indicating that these differences are symptoms of different choices of the generator parameters, which would have been tuned in less extreme regions of phase-space given by previous Z + jets measurements. Measurements in the regions of phase space used as control regions by this and related searches have not been made directly in the past. Therefore, the measurements provided by this paper can be used to help validate new generator versions and parameter choices in this region of phase space, and thus help improve the modelling of the background for future searches.