1 Introduction

Massive coloured particles decaying into quarks and gluons are predicted in several extensions of the Standard Model (SM). At hadron colliders, the search for new phenomena in fully hadronic final states, without missing transverse momentum, is experimentally challenging due to the very large SM multijet production cross-section. This paper describes a search for pair-produced particles each decaying into two jets using 36.7 fb\(^{-1}\) of \(\sqrt{s}=13\) \(\text {TeV}\) proton–proton (pp) collision data recorded in 2015 and 2016 by the ATLAS experiment at the Large Hadron Collider (LHC).

Supersymmetry (SUSY) [1,2,3,4,5,6,7] is a generalisation of the Poincaré symmetry group that relates fermionic and bosonic degrees of freedom. In the generic superpotential, Yukawa couplings can lead to baryon- and lepton-number violation:

$$\begin{aligned} \mathcal {W}_{\text {RPV}} = \lambda _{ijk}L_iL_j\overline{E}_k+\lambda '_{ijk}L_i Q_j\overline{D}_k+\lambda ''_{ijk}\overline{U}_i\overline{D}_j\overline{D}_k + \kappa _iL_iH_u, \end{aligned}$$

where i, j, and k are quark and lepton generation indices. The \(L_i\) and \(Q_i\) represent the lepton and quark \(\mathrm {SU}(2)_\mathrm {L}\) doublet superfields and \(H_u\) the Higgs superfield that couples to up-type quarks. The \(\bar{E}_i\), \(\bar{D}_i\), and \(\bar{U}_i\) are the lepton, down-type quark and up-type quark \(\mathrm {SU}(2)_\mathrm {L}\) singlet superfields, respectively. For each term the couplings are \(\lambda \), \(\lambda '\), \(\lambda ''\), as well as \(\kappa \) which is a dimensional mass parameter. The \(\lambda \) and \(\lambda ''\) couplings are antisymmetric in the exchange of \(i\rightarrow j\) and \(j\rightarrow k\), respectively. While these terms in many scenarios are removed by imposing an additional \(Z_2\) symmetry (R-parity) [8], the possibility that at least some of these R-parity-violating (RPV) couplings are not zero is not ruled out experimentally [9, 10]. This family of models leads to unique collider signatures which can escape conventional searches for R-parity-conserving SUSY.

Naturalness arguments [11, 12] suggest that higgsinos and top squarksFootnote 1 (stops) should be light, with masses below a \(\text {TeV}\) [13, 14]. Third-generation squarks in R-parity-conserving scenarios, and top squarks in particular, have been the subject of a thorough programme of searches at the LHC [15,16,17,18,19,20,21,22].

If the top squark decays through RPV couplings, however, the existing bounds on its mass can be significantly relaxed [23,24,25,26]. Indirect experimental constraints [27] on the sizes of each of the \(\lambda ''\) couplings are primarily valid for low squark mass and for first- and second-generation couplings.

This search targets a model where the top squark is the lightest supersymmetric particle and decays through baryon-number-violating RPV \(\lambda ''\) couplings, \(\tilde{t}\rightarrow \bar{q}_j\bar{q}_k\). The couplings are assumed to be sufficiently large for the decays to be prompt, but small enough to neglect the single-top-squark resonant production through RPV couplings. Top squarks are then produced through strong interactions with cross-sections that do not depend on the specific assumptions in the SUSY model. For two specific choices of couplings, the process considered is schematically depicted in Fig. 1.

In models with extended SUSY, colour-octet states can arise as scalar partners of a Dirac gluino [28,29,30,31]. These scalar gluons (or sgluons) are mostly produced in pairs, and decay into two quarks or two gluons.

Massive colour octet-resonances, generically referred to as colorons (\(\rho \)) [32, 33] are predicted in a wide range of other theories, including axigluon [34, 35] and topcolor [36], in vector-like confinement models [37, 38] and as Kaluza–Klein excitations of the gluons [39, 40]. Colorons can be pair-produced and decay into two jets, a scenario which leads to a four-jet final state.

Constraints on top squarks decaying through \(\lambda ''\) couplings were first set by the ALEPH experiment at LEP [41], excluding at 95% confidence level (CL) masses below 80 \(\text {GeV}\). The CDF experiment at the Tevatron [42], increased these limits to 100 \(\text {GeV}\). Searches for pair-produced resonances in hadronic final states were performed at the LHC at 7 \(\text {TeV}\) and 8 \(\text {TeV}\) of centre-of-mass energy by both the ATLAS [43, 44] and CMS experiments [45, 46]. For decays including heavy-flavour jets in the final state, exclusion limits at 95% CL on the mass of the top squark in the ranges \(100\, \text {GeV}\le m_{\tilde{t}}\le 320\) \(\text {GeV}\) and \(200\, \text {GeV}\le m_{\tilde{t}}\le 385\) \(\text {GeV}\) have been reported by ATLAS [44] and CMS [46], respectively.

Fig. 1
figure 1

Diagrams depicting the direct pair-production of top squarks through strong interactions, with decays into a d- and an s-quark (left) or into a b- and an s-quark (right) through the \(\lambda ''\) R-parity-violating couplings, indicated by the blue dots

2 ATLAS detector

The ATLAS detector [47] is a multi-purpose particle physics detector with a forward-backward symmetric cylindrical geometry and nearly \(4\pi \) coverage in solid angle.Footnote 2 The inner tracking detector consists of pixel and silicon microstrip detectors covering the pseudorapidity region \(|\eta | < 2.5\), surrounded by a transition radiation tracker which provides electron identification in the region \(|\eta | < 2.0\). Starting in Run 2, a new inner pixel layer, the Insertable B-Layer (IBL) [48, 49], has been inserted at a mean sensor radius of 3.3 cm. The inner detector is surrounded by a thin superconducting solenoid providing an 2 T axial magnetic field and by a lead/liquid-argon (LAr) electromagnetic calorimeter covering \(|\eta | < 3.2\). A steel/scintillator-tile calorimeter provides hadronic coverage in the central pseudorapidity range (\(|\eta | < 1.7\)). The endcap and forward regions (\(1.5< |\eta | < 4.9\)) of the hadronic calorimeter are made of LAr active layers with either copper or tungsten as the absorber material. An extensive muon spectrometer with an air-core toroidal magnet system surrounds the calorimeters. Three layers of high-precision tracking chambers provide coverage in the range \(|\eta | < 2.7\), while dedicated fast chambers allow triggering in the region \(|\eta | < 2.4\). The ATLAS trigger system consists of a hardware-based level-1 trigger followed by a software-based high level trigger [50].

3 Data sample

The data used in this analysis were collected by the ATLAS detector in pp collisions at \(\sqrt{s} = 13\) \(\text {TeV}\) at the LHC using a minimum proton bunch crossing interval of 25 ns during 2015 and 2016. In this dataset the mean number of pp interactions per proton bunch crossing is about 23. Events were recorded using a four-jet trigger with transverse momentum (\(p_{\text {T}}\)) thresholds of 100 GeV for each jet at the high-level trigger, which is fully efficient after the analysis selection requirements are applied. After requiring quality criteria for the beam, the data and the detector condition, the available dataset corresponds to an integrated luminosity of 36.7 fb\(^{-1}\) with an uncertainty of \(\pm 2.1\%\) for the 2015 data and \(\pm 3.4\%\) for the 2016 data. The uncertainty in the integrated luminosity is obtained from a calibration of the luminosity scale using a pair of beam-separation scans performed in August 2015 and June 2016, following a methodology similar to that detailed in Ref. [51].

4 Simulated samples

The dominant background from SM multijet production is estimated with a data-driven technique, while Monte Carlo (MC) simulated events are used to estimate the contribution of the \(t\bar{t}\) background, to model the signals and to establish and validate the background estimation method.

The response of the detector was simulated [52] using either a GEANT4 simulation [53] or a fast parameterised simulation [54] of the calorimeter response and GEANT4 for everything else. To account for additional pp interactions in the same and nearby bunch crossings (pile-up), a set of minimum-bias interactions was generated using Pythia 8.186 [55] with the A2 set of parameters (tune) [56] and the MSTW2008LO [57, 58] parton distribution function (PDF) set and was superimposed on the hard scattering events. The EvtGen v1.2.0 program [59] was used to simulate properties of bottom and charm hadron decays for all samples. Corrections were applied to the simulated events to account for differences between data and simulation for the efficiency of identifying jets originating from the fragmentation of b-quarks, together with the probability for mistagging light-flavour and charm-quark jets.

Background samples of multijet production were simulated with \(2\rightarrow 2\) matrix elements (ME) at leading order (LO) using the Pythia 8.186 event generator. The renormalisation and factorisation scales were set to the average \(p_{\text {T}}\) of the two leading jets. The ATLAS A14 tune [60] of parton shower and multiple parton interaction parameters was used together with the NNPDF23LO PDF set [61].

Top-pair production events were simulated using the Powheg-Box v2 [62] generator with the CT10 PDF set. The top mass was set to 172.5 \(\text {GeV}\). The \(h_{\mathrm {damp}}\) parameter, which regulates the transverse momentum of the first extra gluon emission beyond the Born configuration (and thus controls the transverse momentum of the \(t\bar{t}\) system), was set to the mass of the top quark. The parton shower, hadronisation, and underlying event were simulated using Pythia 6.428 [63] with the CTEQ6L1 PDF set and the corresponding Perugia 2012 tune (P2012) [64]. The sample was normalised using the next-to-next-to-leading-order (NNLO) cross-section including the resummation of soft gluon emission at next-to-next-to-logarithmic (NNLL) accuracy using Top++2.0 [65].

The search considers three benchmark signals: the pair production of top squarks, colorons and sgluons with decays into two jets for each resonance.

Signal samples were generated using MG5_aMC@NLO [66] v2.2.3 interfaced to Pythia 8.186 with the A14 tune for the modelling of the parton shower, hadronisation and underlying event. The ME calculation was performed at leading order and, for the top squark signal, includes the emission of up to two additional partons. The merging with the parton shower was done using the CKKW-L [67] prescription, with a merging scale set to one quarter of the pair-produced resonance mass. The PDF set used for the generation is NNPDF23LO. For the top squark signal generation all the non-SM particle masses were set to 5 \(\text {TeV}\) except for the top squark mass (\(m_{\tilde{t}}\)) itself. The top squark was decayed in Pythia 8 assuming a 100% branching ratio into \(\bar{b}\bar{s}\). Its width is expected to be small, and negligible with respect to the detector resolution. This set of samples is also used to interpret the analysis for the case where both top squarks decay into light quarks, since the analysis is not sensitive to the flavour content of the jets. The top squark pair-production cross-sections were calculated at next-to-leading order (NLO) in the strong coupling constant, adding the resummation of soft gluon emission at next-to-leading-logarithmic accuracy [68,69,70]. The nominal cross-section and its uncertainty were taken from an envelope of cross-section predictions using different PDF sets and factorisation and renormalisation scales, as described in Ref. [71]. The coloron samples were generated with the model described in Ref. [72], where the couplings of the vector colour octet to all particles except light quarks were set to zero. The LO cross-sections from the event generator were used. The coloron samples are also used to interpret the result in the context of sgluon pair-production, where they are scaled to the sgluon cross-section computed at NLO with MG5_aMC@NLO [73, 74]. The sgluons are assumed to decay into two gluons, which in this analysis are not distinguished from quark-initiated jets.

5 Event reconstruction

Candidate jets are reconstructed from three-dimensional topological energy clusters [75] in the calorimeter using the anti-\(k_t\) jet algorithm [76], as implemented in the FastJet package [77], with a radius parameter of 0.4. Each topological cluster is calibrated to the electromagnetic energy scale prior to jet reconstruction. The reconstructed jets are then calibrated to the particle level by the application of a jet energy scale (JES) calibration derived from simulation and in situ corrections based on 13 \(\text {TeV}\) data [78,79,80]. The TightBad cleaning quality criteria [81] are imposed to identify jets arising from non-collision sources or detector noise. Any event containing at least one jet failing quality requirements with \(p_{\text {T}} {}>20\) \(\text {GeV}\) is removed.

Jets containing b-hadrons (b-jets) are tagged by a multivariate algorithm (MV2c10) using information about the impact parameters of inner detector tracks associated with the jet, the presence of displaced secondary vertices, and the reconstructed flight paths of b- and c-hadrons inside the jet [82]. A working point with a 77% efficiency, as determined in a simulated sample of \(t\bar{t}\) events, is chosen. The corresponding rejection factors against simulated jets originating from c-quarks and from light quarks or gluons are 4.5 and 130, respectively [83].

6 Event selection

Each event is required to have a reconstructed primary vertex with at least two associated tracks with \(p_{\text {T}} > 400\) \(\text {MeV}\) and a position consistent with the beamspot envelope. If more than one such vertex is found, the vertex with the largest \(\sum p_{\text {T}} ^2\) of the associated tracks is chosen.

The final state under consideration consists of four jets forming two pairs originating from a pair of equal-mass resonances. After the trigger requirement, only events with at least four reconstructed jets with \(p_{\text {T}} > 120\) \(\text {GeV}\) and \(|\eta | < 2.4\) are retained in the analysis.

The analysis strategy exploits the case where the resonances are produced with a significant transverse momentum. As a result the decay products are expected to be close to each other. Taking advantage of this property, candidate resonances are constructed by pairing the four leading jets in the event. Two jet pairs are identified by the following quantity:

$$\begin{aligned} \Delta R_{\text {min}} = \text {min}\left\{ \sum _{i=1,2} |\Delta R_{i} - 1| \right\} , \end{aligned}$$

where \(\Delta R_{i}\) is the angular distance between the two jets for the \(i{\mathrm {th}}\) pair and the sum is over the two pairs of dijets. The offset of \(-1\) is chosen to maximise the signal efficiency for the masses of interest while minimising the effects of soft jets from radiated gluons being recombined with their parent jets in multijet topologies.

The above criteria define the analysis preselection. Additional requirements are applied to further enhance the signal purity. These are based on four discriminating variables established from simulation studies and previous ATLAS searches [38, 43, 44].

To reduce the non-resonant multijet background, for which the pairing efficiency is expected to be poor, a quality criterion is applied to the pairing metric. Resonances with higher masses are produced with a lower boost, and their decay products are less collimated. To compensate for the larger (smaller) angular separation between the jets at high (low) mass this requirement is made dependent on the average reconstructed mass of the two resonance candidates in the event, \(m_{\text {avg}}\). The event is discarded if the best combination of the four leading jets satisfies:

$$\begin{aligned} \Delta R_{\mathrm {min}}&> -0.002 \cdot (m_{\mathrm {avg}}/\text {GeV}-225) + 0.72 \quad \mathrm {if}\; m_{\mathrm {avg}} \le 225~\text {GeV}, \\ \Delta R_{\mathrm {min}}&> +0.0013 \cdot (m_{\mathrm {avg}}/\text {GeV}-225) + 0.72 \quad \mathrm {if}\; m_{\mathrm {avg}} > 225~\text {GeV}. \end{aligned}$$

After boosting the system formed by the two resonances into its centre-of-mass frame, the magnitude of the cosine of the angle that either of them forms with the beamline is denoted as \(|\cos (\theta ^*)|\). Background jets from multijet production frequently originate from t-channel gluon exchange and are preferentially produced in the forward region, with \(|\cos (\theta ^*)|\) close to one. Jets originating from the signal are instead expected to be more central and lead to small \(|\cos (\theta ^*)|\) values.

Since the two reconstructed resonances are expected to have equal mass, their mass difference is a powerful discriminant between signal and background. The mass asymmetry (\(\mathcal {A}\)) is defined as:

$$\begin{aligned} \mathcal {A} = \frac{|m_{1} - m_{2}|}{m_{1} + m_{2}}, \end{aligned}$$

where \(m_{1}\) and \(m_{2}\) are the invariant masses of the two reconstructed dijet pairs. The value of \(\mathcal {A}\) is expected to peak at zero for well-paired signal events and to have larger values for background events

Fig. 2
figure 2

The distributions of the (a) smallest angular separation between the two jets in a pair (\(\Delta R_{\text {min}}\)), the (b) mass asymmetry (\(\mathcal {A}\)), the (c) pair production angle \(|\cos (\theta ^*)|\) and the (d) multiplicity of b-tagged jets. The observed data (black dots) are compared with the distributions expected from a top squark with a mass of 250 \(\text {GeV}\) (solid blue line) or 500 \(\text {GeV}\) (azure dotted line) and a coloron with a mass of 1500 \(\text {GeV}\) (red dashed line). The distributions are normalised to unity and shown at preselection, after the requirement of four jets paired into two candidate resonances

The distributions of \(\Delta R_{\mathrm {min}}\), \(\mathcal {A}\) and \(|\cos (\theta ^*)|\) after preselection are shown for data, a top squark sample with a mass of \(m_{\tilde{t}}=500\) \(\text {GeV}\) and a coloron sample with mass \(m_{\rho }=1500\) \(\text {GeV}\) in Fig. 2a–c. Because of the very small expected signal purity (below 2%) before additional selection criteria are applied the data distributions can be viewed as representative of the expected background. Two additional requirements, \(\mathcal {A} < 0.05\) and \(|\cos (\theta ^*)|<0.3\), define the inclusive signal region (SR) selection, targeting resonance decays into light quark or gluon jets. The selections are determined in an optimisation procedure that maximises the expected signal significance.

When the dominant RPV couplings involve third-generation quarks (\(\lambda ^{''}_{3i3}\)), a b-quark is expected from each of the top squark decays. A dedicated b-tagged SR selection is used for this scenario. On top of the requirements applied in the inclusive selection it requires at least two b-tagged jets to be present in the event, which significantly reduces the multijet background. The distribution of the number of b-tagged jets after pairing the four jets into candidate resonances is shown for data and two top squark signals with masses of 250 and 500 \(\text {GeV}\) in Fig. 2d. An additional factor of about two in background reduction is gained by requiring each of the two b-jets to be associated with a different reconstructed resonance. This is particularly effective in reducing the contribution of \(g\rightarrow b\bar{b}\) splittings, where the two b-jets are typically very collimated.

The final analysis discriminant is the average mass of the two reconstructed resonances:

$$\begin{aligned} m_{\mathrm {avg}}=\frac{1}{2}(m_1+m_2). \end{aligned}$$

A peak in \(m_{\text {avg}}\) at a mass of about that of the resonance is expected for the signal, over a non-peaking background from multijet processes. Figure 3 shows the expected \(m_{\text {avg}}\) distribution for signal samples with different masses. For each mass hypothesis a counting experiment is performed in a window of the \(m_{\text {avg}}\) variable optimised to maximise the expected signal significance. The windows range from a 10 \(\text {GeV}\) width for a 100 \(\text {GeV}\) top squark to a 200 \(\text {GeV}\) width for a 1500 \(\text {GeV}\) coloron. The mass window for the highest target mass considered, of 2000 \(\text {GeV}\), has no upper edge. When the mass difference between two signal samples is smaller than the experimental resolution, their selected mass windows partially overlap.

Fig. 3
figure 3

Distribution of the average mass, \(m_{\text {avg}}\), in the inclusive signal region for simulated top squark signals with \(m_{\tilde{t}}=250\), 500, and 750 and coloron signals with \(m_{\rho }=1000\), 1250, and 1500 \(\text {GeV}\)

Fig. 4
figure 4

The acceptance times efficiency (Acc.\(\times \epsilon \)) of the a inclusive and b b-tagged signal region selection as a function of the resonance mass, m, before and after the mass window requirements are applied. Top squark signals are indicated by the blue triangles, coloron by the red squares. The statistical uncertainties are indicated by vertical bars

Table 1 MC predictions of the number of signal events corresponding to 36.7 fb\(^{-1}\) of data after applying each of the event selection requirements, except for the mass window. Top squark masses of \(m_{\tilde{t}}=100\) \(\text {GeV}\) and \(m_{\tilde{t}}=500\) \(\text {GeV}\), and a coloron mass of 1500 \(\text {GeV}\) are shown. The statistical uncertainty of the MC simulation is shown for each selection

The MC predictions of signal event yields in 36.7 fb\(^{-1}\) of data are shown in Table 1 after each different requirement of the event selection is applied. The acceptance times efficiency of the inclusive and b-tagged signal region selections as a function of the signal mass are shown before and after applying the \(m_{\text {avg}}\) mass window requirement in Fig. 4. The acceptance of the signal region selections increases for large masses due to the four jets from the signal having a larger \(p_{\text {T}}\). However, as the jet pairing does not always correctly assign the resonance candidates for high masses, the signal has a tail extending to low \(m_{\text {avg}}\) values, degrading the efficiency of the mass window selection.

7 Background estimation

The dominant background from multijet production is estimated directly from data, with a method that predicts both the normalisation and the shape of the \(m_{\text {avg}}\) distribution. In the b-tagged selection, for \(m_{\text {avg}}\) below 200 \(\text {GeV}\), the \(t\bar{t}\) contribution becomes significant, and is estimated from simulation.

Fig. 5
figure 5

Definition of the control and validation regions in the \(\mathcal {A}\) and \(|\cos (\theta ^*)|\) plane used to estimate the multijet background

For the inclusive selection, the \(m_{\text {avg}}\) distribution for the background is obtained from data. For each \(m_{\text {avg}}\) bin the data sample is divided into four regions: one region where the signal region selection is applied (D) and three background-dominated control regions (A, C and F). The variables used to define the different regions, summarised in Fig. 5, are \(\mathcal {A}\) and \(|\cos (\theta ^*)|\). Provided the two variables defining the regions are uncorrelated, and signal leakage in the background-dominated regions can be neglected, the amount of background in the region of interest D can be predicted from the observed numbers of events in the control regions as \(N_\mathrm {D} = N_\mathrm {A} \times N_\mathrm {F}/N_\mathrm {C}\). The linear correlation between the \(|\cos (\theta ^*)|\) and \(\mathcal {A}\) variables is evaluated in data and simulated multijet samples, where it amounts to 1.8 and 2.2%, respectively. Significant correlations are observed in data at large \(m_{\text {avg}}\) and high \(\mathcal {A}\) values; to reduce their impact on the background estimate the \(\mathcal {A}\)\(|\cos (\theta ^*)|\) plane is restricted to \(0.0<\mathcal {A}<0.7\) and \(0.0<|\cos (\theta ^*)|<0.7\). Two additional regions (B and E) are defined in the \(\mathcal {A}\)\(|\cos (\theta ^*)|\) plane. The validation region (VR), region E, is used to test the performance of the data-driven method and assign an uncertainty to the background estimate. The validation region is defined with the same selections as for the signal region, but with the asymmetry requirement changed from \(\mathcal {A} < 0.05\) to \(0.05< \mathcal {A} < 0.15\). The background contribution in the VR is estimated by \(N_\mathrm {E} = N_\mathrm {B} \times N_\mathrm {F}/N_\mathrm {C}\). In the inclusive selection the data-driven estimate also accounts for the contribution from the \(t\bar{t}\) production, which amounts to less than 1% of the total background for \(m_{\mathrm {avg}}<200\) \(\text {GeV}\), and is negligible above.

For the b-tagged selection, where the background is relatively small, the signal contamination in region \(\mathrm {A}\) can be significant and potentially bias the result of the background estimate. The multijet background for this selection is thus estimated in two steps. The shape of the \(m_{\text {avg}}\) distribution is first predicted in a region with a b-tag veto (zero-tag) and then extrapolated to the b-tagged signal region. The \(m_{\text {avg}}\) distribution in the zero-tag region is obtained with a data-driven estimate, analogously to the inclusive selection. The zero-tag prediction is then extrapolated to the b-tagged selection by means of projection factors computed bin by bin in \(m_{\text {avg}}\), similarly to the approach described in Ref. [44]. The projection factors, for a given \(m_{\text {avg}}\) bin and region in the \(\mathcal {A}\)\(|\cos (\theta ^*)|\) plane, are defined as the ratio of the numbers of events with two b-tags and zero b-tags, \(N_{\mathrm {two-}b\mathrm {-tags}}/N_{\mathrm {zero-}b\mathrm {-tags}}\), within that region. The method assumes the projection factors to be constant across the \(\mathcal {A}\)\(|\cos (\theta ^*)|\) plane. They are evaluated in region F, where a negligible signal contamination is expected. The contributions from multi-jet and \(t\bar{t}\) production scale differently between the zero- and b-tagged selection. Hence, simulated samples are used to subtract the \(t\bar{t}\) contribution in all control regions. The \(t\bar{t}\) estimate in the signal region is then obtained directly from the simulation, considering all relevant modelling and experimental uncertainties.

Table 2 Observed numbers of events and the predicted \(t\bar{t}\) contributions in each of the regions used in the background estimate, for each b-tag multiplicity. The expected fractional signal contributions are shown for the mass windows corresponding to \(m_{\tilde{t}}=125\), 250, 500, and 800 \(\text {GeV}\). For the \(m_{\tilde{t}}=125\) and 250 \(\text {GeV}\) mass windows the fractions of \(t\bar{t}\) are also shown. The \(t\bar{t}\) systematic uncertainties include both the detector-level uncertainties and the theoretical uncertainties, as described in Sect. 8
Fig. 6
figure 6

The \(m_{\text {avg}}\) spectrum in the inclusive (left) and b-tagged (right) validation regions. The data (black points) are compared to the total background prediction (red line) estimated with the data-driven method. The fraction of background coming from top-pair production is shown in orange. The statistical uncertainties of the prediction are shown by the grey hatched band

The observed number of events in each of the regions used in the background estimate before the mass window requirements are applied, together with the expected signal contamination in few representative mass windows, are shown for both the inclusive and b-tagged selections in Table 2. The \(m_{\text {avg}}\) distribution in the validation region for the inclusive and b-tagged signal regions is shown in Fig. 6. Within the statistical uncertainties the method reproduces both the normalisation and the shape of \(m_{\text {avg}}\) in the VRs. The level of agreement observed in the VRs is used to derive a systematic uncertainty in the background estimate in the SR. In each \(m_{\text {avg}}\) mass window the difference between the observed data and the estimation in the VR (non-closure) is computed. The larger of the observed non-closure in the VR and the statistical uncertainty of the data-driven method is assigned as an uncertainty in the background estimates. To reduce the effect of statistical fluctuations in the non-closure and avoid quoting an unphysically small value of the systematic uncertainty for the mass windows where it changes sign, this uncertainty is further smoothed as a function of \(m_{\text {avg}}\). The Nadaraya–Watson [84, 85] kernel regression estimate is used for the smoothing, with a bandwidth of 500 \(\text {GeV}\) (meaning that the quartiles of the kernels are placed at \(\pm 125\) \(\text {GeV}\)). The uncertainties assigned to the background estimate in the inclusive and b-tagged signal regions are summarised in Fig. 7.

Fig. 7
figure 7

The uncertainty in the data-driven background estimate in the inclusive (left) and b-tagged (right) signal regions, computed in the \(m_{\text {avg}}\) mass windows defined for different target masses. The uncertainty arising from the non-closure in the validation region is shown with a red short-dashed line and compared with the statistical uncertainty of the data-driven prediction shown as an orange dashed line. The additional uncertainty in the MC estimate of the top background in the b-tagged signal region is shown as a dotted blue line. The total uncertainty, obtained by adding in quadrature the different components, is indicated by a black solid line

8 Systematic uncertainties

While the multijet background uncertainties pertain primarily to the estimation method itself, the top background and the signals are also affected by uncertainties related to the description of detector effects and to the physics modelling of the MC simulation.

The dominant detector-related systematic effects are due to the uncertainties in the jet energy scale [80] and resolution [86] and from the b-tagging efficiency and mistag rate [83].

Since MC simulation is used to determine the contribution from top events in the b-tagged signal region, systematic uncertainties related to the choice of MC generator for the process need to be estimated. These are evaluated by comparing the nominal samples to additional samples with systematic variations. A modelling uncertainty is derived by comparing the predictions of the nominal sample with a sample produced with Powheg interfaced with Herwig++ 2.7.1, or with MG5_aMC@NLO interfaced with Herwig++. In addition, the difference in the prediction obtained by modifying the parton-shower intensity and the \(h_{\mathrm {damp}}\) parameter in the nominal sample is taken as an uncertainty. The \(t\bar{t}\) systematic uncertainty on the total background is as large as 20% for \(m_{\text {avg}}\) below 200 \(\text {GeV}\), becoming negligible above 200 \(\text {GeV}\).

The total detector-related uncertainties in the signal are about 10% for the inclusive SR and about 15% for the two-b-tagged SR. For top squark production the nominal signal cross-sections and their uncertainties are taken from an envelope of cross-section predictions derived using different PDF sets and different factorisation and renormalisation scales, as described in Ref. [71]. The theoretical uncertainties in the acceptance of the signal simulation include variations of the renormalisation and factorisation scales, the CKKW-L merging scales, and the value of the strong coupling constant in MG5_aMC@NLO as well as parton shower uncertainties in Pythia 8 evaluated from variations of the A14 parameter set. After normalising the samples using the same cross-section, the difference between the yields from the nominal and varied samples in the mass window, which is typically below 1%, is considered as an uncertainty.

9 Results and interpretation

The \(m_{\text {avg}}\) distributions in the inclusive and b-tagged signal regions are shown in Fig. 8. Agreement is observed between data and the expected background. The expected numbers of background and signal events in the SR and their uncertainties are reported for the mass windows defined for top squark and coloron signals in Tables 3 and 4, respectively, together with the observed events in data. Table 5 presents the numbers in the top squark mass windows of the two-b-tagged signal region.

To estimate the compatibility of the data with a generic resonance mass hypothesis, the \(m_{\text {avg}}\) distribution is scanned in 12.5 \(\text {GeV}\) steps. The \(m_{\text {avg}}\) window for an arbitrary mass is obtained from a linear fit to the lower and upper edges of the windows obtained for the simulated signal masses. For each mass a background \(p_0\)-value is computed for the inclusive and b-tagged signal regions. The largest deviation is found in the b-tagged signal region for a mass of 463 \(\text {GeV}\), corresponding to a local \(p_0\) value of 0.05.

The expected \(p_0\) values in each mass window are also evaluated for potential signals. At least three-standard-deviation \((3\sigma )\) signal sensitivity is expected for top squark masses up to 350 \(\text {GeV}\) with the inclusive signal region and 450 \(\text {GeV}\) with the b-tagged signal region. For colorons a greater than \(3\sigma \) sensitivity is expected for masses up to 1400 \(\text {GeV}\).

Table 3 Observed numbers of events in the data, \(N_{\mathrm {Data}}\), the estimated numbers of background events, \(N_{\mathrm {Bkg}}\), and the expected numbers of top squark signal events, \(N_{\mathrm {Sig}}\), in the top squark mass windows of the inclusive signal region. Separate statistical and systematic uncertainties are given
Table 4 Observed numbers of events in the data, \(N_{\mathrm {Data}}\), the estimated numbers of background events, \(N_{\mathrm {Bkg}}\), and the expected numbers of coloron signal events, \(N_{\mathrm {Sig}}\), in the coloron mass windows of the inclusive signal region. Separate statistical and systematic uncertainties are given
Table 5 Observed numbers of events in the data, \(N_{\mathrm {Data}}\), the estimated numbers of background events, \(N_{\mathrm {Bkg}}\), and the expected numbers of top squark signal events, \(N_{\mathrm {Sig}}\), in the top squark mass windows of the b-tagged signal region. Separate statistical and systematic uncertainties are given
Fig. 8
figure 8

The \(m_{\text {avg}}\) spectrum in the inclusive (left) and b-tagged (right) signal regions. The data (black points) are compared to the total background prediction (red line) estimated with the data-driven method. The fraction of background coming from top-pair production is shown in orange. The statistical uncertainties of the prediction are shown by the grey hatched band. Signals of different masses are overlaid in different colours

In the absence of a statistically significant excess in data, exclusion limits are derived for the investigated signal models. The inclusive signal region is used to set a limit on top squark, sgluon and coloron production with decays into a pair of jets, while the b-tagged signal region is used to interpret the search for top squark pair production with decays into a b- and a light-quark jet. A profile likelihood ratio combining Poisson probabilities for signal and background is computed to determine the 95% CL interval for compatibility of the data with the signal-plus-background hypothesis (\(\mathrm {CL}_\mathrm {s+b}\)) [87]. A similar calculation is performed for the background-only hypothesis (\(\mathrm {CL}_\mathrm {b}\)). From the ratio of these two quantities, the confidence level for the presence of a signal (\(\mathrm {CL}_\mathrm {s}\)) is determined [88]. Systematic uncertainties are treated as nuisance parameters and are assumed to follow Gaussian distributions. The results are evaluated using pseudo-experiments. This procedure is implemented using a software framework for statistical data analysis, HistFitter [89]. The observed and expected 95% CL upper limits on the allowed cross-sections are shown in Fig. 9. For top squark decays into two quarks, the expected and observed mass range exclusions are between 100 and 430 \(\text {GeV}\) and between 100 and 410 \(\text {GeV}\), respectively. This exclusion does also apply to the pair-production of other squarks, decaying, for example, to a d- and a u-quark. If the top squark decay is into a b-quark and a light-quark, masses between 100 and 530 \(\text {GeV}\) are expected to be excluded, with the observed exclusion ranging from 100 to 470 \(\text {GeV}\) and from 480 to 610 \(\text {GeV}\). Below top squark masses of about 200 \(\text {GeV}\) the signal acceptance rapidly drops due to the trigger and jet requirements, and the analysis sensitivity does not surpass the 8 \(\text {TeV}\) result, which was specifically optimised for low-mass signals. Pair-produced scalar gluons with decays into two gluons are excluded up to a mass of 800 \(\text {GeV}\). Pair-produced colorons coupling only to light quarks are excluded up to a mass of 1500 \(\text {GeV}\).

Fig. 9
figure 9

The 95% CL upper limit on the \(\sigma \times B\) value compared to the theoretical cross-section for the direct pair-production of top squarks with decays into a \(\bar{q}\bar{q'}\) or b \(\bar{b}\bar{s}\) and c high-mass colorons decaying into qq and sgluons decaying into gg. The dashed black and solid red lines show the 95% CL expected and observed limits respectively, including all uncertainties except the theoretical signal cross-section uncertainty. The solid green (yellow) band around the expected limit shows the associated \(\pm 1\sigma \) (\(\pm 2\sigma \)) ranges. The shaded coloured cross-section band indicates the \(\pm 1\sigma \) variations due to theoretical uncertainties in the signal production cross-section given by renormalisation and factorisation scale and PDF uncertainties. The region of low top squark mass not shown in the plot is excluded by Refs. [41, 42]

10 Conclusion

A search is presented for the pair production of coloured resonances, each decaying into two jets. The analysis uses 36.7 fb\(^{-1}\) of \(\sqrt{s}\) = 13 TeV pp collision data recorded by the ATLAS experiment at the LHC in 2015 and 2016. An inclusive selection and a selection with two b-tagged jets in the event are defined, and counting experiments are performed in windows of the average mass of the two resonance candidates. No significant deviation from the background prediction is observed. The results are interpreted in a SUSY simplified model with a top squark as the lightest supersymmetric particle, which is pair-produced and decays promptly into two quarks through R-parity-violating couplings. For decays into two quarks, top squark masses in the range \(100\,\text {GeV}<m_{\tilde{t}}<410\,\text {GeV}\) are excluded at 95% CL. If the top squark decays into a b-quark and a light quark, masses in the ranges \(100\,\text {GeV}<m_{\tilde{t}}<470\,\text {GeV}\) and \(480\,\text {GeV}<m_{\tilde{t}}<610\,\text {GeV}\) are excluded at 95% CL. Limits on the pair production of scalar gluons with decays into two gluons reach masses of 800 \(\text {GeV}\). Vector colour-octet resonances coupling only to light quarks are excluded up to masses of 1500 \(\text {GeV}\). The results improve upon previous Run 1 searches and extend the constraints on top squark masses.