1 Introduction

Pairs of highly energetic jets (dijets) are produced at high rates in proton–proton collisions at the CERN LHC through pointlike scattering of quarks and gluons. Despite its enormous success, the shortcomings of the standard model (SM) are well known. Many theories of physics beyond the standard model (BSM) that alter the interaction of quarks and gluons from that predicted by perturbative quantum chromodynamics (QCD) give rise to narrow or wide resonances or even to nonresonant dijet signatures. Examples that have received widespread attention include models with dark matter (DM) [1,2,3,4,5], quark compositeness [6,7,8], extra spatial dimensions [9, 10], and quantum black holes [11,12,13,14,15]. Resonances with an intrinsic width of the order of the experimental resolution can be constrained by searches in the dijet invariant mass spectrum [16,17,18]. These searches, however, are not very sensitive to wide resonances or nonresonant signatures; a more effective strategy to constrain such signatures is the study of dijet angular distributions [19].

The angular distribution of dijets relative to the beam direction is sensitive to the dynamics of the scattering process. Furthermore, since the angular distributions of the dominant underlying QCD processes of \(\mathrm {q}\mathrm {g}\rightarrow \mathrm {q}\mathrm {g}\), \(\mathrm {q}\mathrm {q}' \rightarrow \mathrm {q}\mathrm {q}'\), \(\mathrm {q}\mathrm {q}\rightarrow \mathrm {q}\mathrm {q}\), \(\mathrm {g}\mathrm {g}\rightarrow \mathrm {g}\mathrm {g}\), are all similar [20], the dijet angular distribution is insensitive to uncertainties in the parton distribution functions (PDFs). The dijet angular distribution is typically expressed in terms of \(\chi _{\text {dijet}} = \exp (|y_1 -y_2 |)\), where \(y_1\) and \(y_2\) are the rapidities of the two jets with the highest transverse momentum \(p_{\mathrm {T}}\) (the leading jets). For collinear massless parton scattering, \(\chi _{\text {dijet}}\) takes the form \(\chi _{\text {dijet}} =(1+|\cos \theta ^{*} |)/(1-|\cos \theta ^{*} |)\), where \(\theta ^{*}\) is the polar scattering angle in the parton-parton center-of-mass (CM) frame. The choice of \(\chi _{\text {dijet}}\), rather than \(\theta ^{*}\), to measure the dijet angular distribution is motivated by the fact that in Rutherford scattering, where only t-channel scattering contributes to the partonic cross section, the \(\chi _{\text {dijet}}\) distribution is independent of \(|y_1 -y_2 |\) [20]. In contrast, BSM processes may have scattering angle distributions that are closer to being isotropic than those given by QCD processes and can be identified by an excess of events at small values of \(\chi _{\text {dijet}}\). Previous measurements of dijet angular distributions at the LHC have been reported by the ATLAS [17, 21,22,23,24,25] and CMS [26,27,28,29] Collaborations.

In a simplified model of interactions between DM particles and quarks  [1,2,3,4, 30, 31], the spin-1 (vector or axial-vector) DM mediator particle with unknown mass \(M_{\text {Med}}\) is assumed to decay only to pairs of quarks or pairs of DM particles, with mass \(m_{\mathrm {DM}}\), and with a universal quark coupling \(g_{\mathrm {\mathrm {q}}}\) and a DM coupling \(g_{\mathrm {DM}}\). In this model, the relative width of the DM mediator increases monotonically with increasing \(g_{\mathrm {\mathrm {q}}}\). In a scenario where \(g_{\mathrm {\mathrm {q}}} =0.25\) and in which the relative widths for vector and axial-vector mediators in the dark matter decay channels are negligible, values of \(M_{\text {Med}}\) below 3.0TeVwere excluded by narrow dijet resonance searches [17, 18]. A search for narrow and broad dijet resonances set constraints on mediator widths up to 30% (\(g_{\mathrm {\mathrm {q}}} <0.75\)) and masses up to 4 TeV [32]. Searches for invisible particles produced in association with quarks or bosons [33,34,35] have excluded vector and axial-vector mediators below 1.8 (2.1)TeVfor \(g_{\mathrm {\mathrm {q}}} =0.25\) (\(g_{\mathrm {\mathrm {q}}} =1.0\)) and \(g_{\mathrm {DM}} =1.0\) [34].

A common signature of quark compositeness [6,7,8], at energies well below the characteristic mass scale \(\varLambda \) for new interactions between quark constituents, is the four-fermion contact interaction (CI). The most stringent limits on quark CIs come from searches in dijet angular distributions at large dijet invariant masses (\(M_{\mathrm {jj}}\)) [17, 29], and in inclusive jet \(p_{\mathrm {T}}\) distributions [36]. The publication from the ATLAS Collaboration [17] provides lower limits on the quark CI scales from 13.1 to 29.5TeV, depending on the details of the model.

The Arkani–Hamed–Dimopoulos–Dvali (ADD) model [9, 10] of compactified large extra dimensions (EDs) provides a possible solution to the hierarchy problem of the standard model. It predicts signatures of virtual graviton exchange that result in a nonresonant enhancement of dijet production in proton–proton collisions, whose angular distribution differs from the predictions of QCD. Signatures from virtual graviton exchange have previously been sought at the LHC in various final states, where the most stringent limits arise from the CMS search with dijet angular distributions [29], which excludes the ultraviolet cutoff in the ADD framework up to 7.9–11.2TeV, depending on the parameterization of the model.

In models with large EDs, the fundamental Planck scale (\({M_\mathrm {Pl}}\)) is assumed to be closer to the electroweak (EW) scale, thereby allowing black hole production at the LHC [11,12,13,14,15]. Semiclassical black holes, which have mass much larger than \({M_\mathrm {Pl}}\), decay into multiple jets through Hawking radiation [37]. Quantum black holes (QBHs), which are produced with mass close to \({M_\mathrm {Pl}}\), decay predominantly into dijets and can be studied using dijet angular distributions [38,39,40]. Recent searches for QBHs with dijet final states at the LHC reported in Refs. [17, 29] exclude QBHs with masses below 8.9TeV.

In this paper, we present a search for new physics, specifically DM mediators, CIs, EDs, and QBHs, using measurements of dijet angular distributions. The signature of the signals can be categorized into nonresonant excesses at high \(M_{\mathrm {jj}}\) as predicted by the CI and ADD models and resonances from the decay of QBHs and DM mediators that could appear across the whole range of the \(M_{\mathrm {jj}}\) spectrum. The searches are performed by comparing detector-level dijet angular distributions with BSM predictions that have been adjusted to include detector resolution effects. This eliminates some systematic uncertainties that are introduced when correcting the dijet angular distributions for detector effects and simplifies the statistical evaluation. The dijet angular distributions are also corrected to particle level to facilitate comparisons with other theoretical predictions and published in HEPData.

2 The CMS detector

The CMS apparatus is based on a superconducting solenoid of 6\(\text {\,m}\) internal diameter, providing an axial field of 3.8\(\text {\,T}\). Within the solenoid and nearest to the interaction point are the silicon pixel and strip trackers. Surrounding the tracker volume are the lead tungstate crystal electromagnetic calorimeter and the brass and scintillator hadron calorimeter. The trackers cover a pseudorapidity region of \(|\eta | < 2.5\) while the calorimeters cover \(|\eta | < 3.0\). In addition, CMS has extensive forward calorimetry, which extends the coverage to \(|\eta |<5.0\). Finally, muons are measured in gas-ionization detectors embedded in the steel flux-return yoke of the solenoid, with a coverage of \(|\eta | < 2.4\). A two-tiered system, with a level-1 trigger followed by a high-level trigger (HLT), is used by CMS to record events of interest [41] for the offline analysis. A more detailed description of the CMS detector, together with a definition of the coordinate system used and the relevant kinematic variables, can be found in Ref. [42].

3 Event selection and data unfolding

Events are reconstructed using a particle-flow algorithm [43] to identify and reconstruct individual particles from each collision by combining information from all CMS subdetectors. Identified particles include charged and neutral hadrons, electrons, muons, and photons. The particles are clustered into jets using the anti-\(k_{\mathrm {T}}\) algorithm [44, 45] with a distance parameter of 0.4. In order to mitigate the effect of additional proton–proton interactions within the same or nearby bunch crossings (pileup) on the jet momentum measurement, the charged hadron subtraction technique [43] is used. Spurious jets from noise or non-collision backgrounds are rejected by applying jet identification criteria [46]. The jet energies are corrected for nonlinear and nonuniform response of the calorimeters through corrections obtained from data and Monte Carlo (MC) simulations [47]. To compare data with theoretical predictions, the same jet clustering algorithm is applied to the generated stable particles (lifetime \(c\tau > 1\,\text {cm} \)) from MC simulations with leading order (LO) pythia 8.212 [48, 49] predictions, and to the outgoing partons from next-to-leading (NLO) predictions.

The events used in this analysis are selected with triggers based upon either jet \(p_{\mathrm {T}}\) or \(H_{\mathrm {T}}\), as measured by the HLT, where \(H_{\mathrm {T}}\) is the scalar sum of the \(p_{\mathrm {T}}\) values of all the jets with \(|\eta | < 3.0\) and \(p_{\mathrm {T}}\) greater than 30GeV. The HLT selection requires having a jet with \(p_{\mathrm {T}} >450\hbox {GeV}\) or an \(H_{\mathrm {T}}\) value of at least 900GeV. The data sample was collected with the CMS detector in 2016 and corresponds to an integrated luminosity of 35.9\(\,\text {fb}^{-1}\)  [50].

In the subsequent offline analysis, events with a reconstructed primary vertex that lies within ± 24 \(\text {\,cm}\) of the detector center along the beam line, and within 2\(\text {\,cm}\) of the detector center in the plane transverse to the beam, are selected. The primary vertex is defined as the reconstructed vertex with the highest sum of the squares of all associated physics objects \(p_{\mathrm {T}}\). The physics objects are the jets returned by the application of the anti-\(k_{\mathrm {T}}\) algorithm to all tracks associated with the vertex, plus the corresponding associated missing transverse momentum, taken as the negative vector sum of the \(p_{\mathrm {T}}\) of those jets.

The two leading jets are used to measure the dijet angular distributions in seven regions of the dijet invariant mass \(M_{\mathrm {jj}}\). The \(M_{\mathrm {jj}}\) regions, in units of TeV, are chosen to be 2.4–3.0, 3.0–3.6, 3.6–4.2, 4.2–4.8, 4.8–5.4, 5.4–6.0, and > 6.0. The highest \(M_{\mathrm {jj}}\) range was chosen to maximize the expected sensitivity to the BSM signals considered. The phase space for this analysis is restricted by the requirements \(\chi _{\text {dijet}} <16\) and \(|y_{\text {boost}} |<1.11\), where \(y_{\text {boost}} = (y_1 + y_2)/2\). This selection and the \(M_{\mathrm {jj}}\) range definition restrict the absolute rapidities \(|y_1 |\) and \(|y_2 |\) of the two highest \(p_{\mathrm {T}}\) jets to be less than 2.5 and their \(p_{\mathrm {T}}\) to be larger than 200GeV. The trigger efficiency for events that satisfy the subsequent selection criteria exceeds 99% in all the \(M_{\mathrm {jj}}\) ranges for the analysis. The observed numbers of events in the analysis phase space for each of the mass ranges are 353025, 71832, 16712, 4287, 1153, 330, and 95. The highest value of \(M_{\mathrm {jj}}\) observed among these events is 8.2TeV.

In this paper, we present dijet angular distributions normalized to unity in each \(M_{\mathrm {jj}}\) range, denoted \((1/\sigma _\text {dijet})(\mathrm {d}\sigma _\text {dijet}/\mathrm {d}\chi _\text {dijet})\), where \(\sigma _\text {dijet}\) is the cross section in the analysis phase space.

Fluctuations in jet response from the resolution in jet \(p_{\mathrm {T}}\) of the detector can cause lower energy jets to be misidentified as leading jets and also result in bin-to-bin event migrations in both \(\chi _{\text {dijet}}\) and dijet mass. The corrections for these effects are obtained from a two-dimensional response matrix that maps the generator-level \(M_{\mathrm {jj}}\) and \(\chi _{\text {dijet}}\) distributions onto the measured values. This matrix is obtained using particle-level jets from the pythia MC event generator that are smeared in \(p_{\mathrm {T}}\) using a double-sided Crystal Ball parameterization [51] of the response. This parameterization takes into account the full jet energy resolution, including non-Gaussian tails, and is derived from the full detector simulation. The width of the Gaussian core in the parameterization is adjusted to account for the difference in resolution observed between data and simulation [47]. The reason for deriving the response matrix from smeared generator-level MC rather than from full detector simulation is that significantly smaller statistical uncertainties can be achieved using the faster code. The measured distributions are unfolded to particle level by inverting the response matrix without regularization, using the RooUnfold package [52]. The unfolding changes the shape of the \(\chi _{\text {dijet}}\) distributions by < 1% and < 8% across \(\chi _{\text {dijet}}\) in the lowest and highest \(M_{\mathrm {jj}}\) ranges, respectively. The fractions of event migrations between mass bins are 15–20% in the lowest \(M_{\mathrm {jj}}\) range and 25–40% in the highest \(M_{\mathrm {jj}}\) range, depending on \(\chi _{\text {dijet}}\) values. The unfolding procedure was tested by splitting the simulation data into independent training and testing samples. The training sample was used to derive a response matrix and the smeared \(\chi _{\text {dijet}}\) distributions from the test sample were unfolded using this response matrix. No significant difference was observed between the generated and unfolded \(\chi _{\text {dijet}}\) distributions in the test sample. The effects of migrations between \(\chi _{\text {dijet}}\) bins are negligible. The unfolding procedure is based on matrix inversion, while the procedure used in previous publications of dijet angular distributions [28, 29] was based on the D’Agostini iterative method [53]. We have compared these two methods by deriving limits from unfolded data, and the limits vary by less than 5%.

4 Theoretical predictions

We compare the unfolded normalized dijet angular distributions with the predictions of perturbative QCD at NLO, available in nlojet++ 4.1.3 [54] in the fastnlo 2.1 framework [55]. EW corrections for dijet production [56] change the predicted normalized distributions by up to 1% (5%) for the lowest \(\chi _{\text {dijet}}\) bins in small (large) values of \(M_{\mathrm {jj}}\). The factorization (\(\mu _{\mathrm {f}}\)) and renormalization (\(\mu _{\mathrm {r}}\)) scales are set to the average \(p_{\mathrm {T}}\) of the two jets, \(\langle p_\mathrm {T} \rangle =(p_{\mathrm {T}} {}_{1}+p_{\mathrm {T}} {}_{2})/2\), and the PDFs are taken from the CT14 set [57]. The use of a more flexible statistical combination of multiple PDF sets as in PDF4LHC15_100 [57,58,59,60,61,62] exhibited small differences as compared to the CT14 PDF set. We evaluated the impact of nonperturbative effects from hadronization and multiple parton interactions on the QCD predictions using pythia with the CUETP8M1 tune [63] and herwig++ 2.7.1 [64] with tune EE5C [65]. The effects are found to be less than 1% and negligible for both MC generators.

The production and decay of the DM mediators in the simplified DM model are generated at LO using MadDM version 2.0.6 [66, 67] at fixed \(g_{\mathrm {DM}}\) and \(m_{\mathrm {DM}}\) values, where \(g_{\mathrm {DM}} = 1.0\) and \(m_{\mathrm {DM}} = 1\hbox {GeV}\). For these values of \(g_{\mathrm {DM}}\) and \(m_{\mathrm {DM}}\), the differences between vector and axial-vector mediators in the cross sections and in the acceptances are negligible in the analysis phase space.

BSM physics signatures from CIs with flavor-diagonal color-singlet couplings among quarks are described by the effective Lagrangian [7, 8]:

$$\begin{aligned} \begin{aligned} \mathcal {L}_{\mathrm {q}\mathrm {q}}&=\frac{2\pi }{\varLambda ^2} [ \eta _{\mathrm {LL}} ({\overline{\mathrm {q}}}_{\mathrm {L}}\gamma ^{\mu }\mathrm {q}_{\mathrm {L}})(\overline{\mathrm {q}} _{\mathrm {L}}\gamma _{\mu }\mathrm {q}_{\mathrm {L}}) \\&\quad +\,\eta _{\mathrm {RR}} (\overline{\mathrm {q}} _{\mathrm {R}}\gamma ^{\mu }\mathrm {q}_{\mathrm {R}})(\overline{\mathrm {q}} _{\mathrm {R}}\gamma _{\mu }\mathrm {q}_{\mathrm {R}}) \\&\quad +\,2\eta _{\mathrm {RL}} (\overline{\mathrm {q}} _{\mathrm {R}}\gamma ^{\mu }\mathrm {q}_{\mathrm {R}})(\overline{\mathrm {q}} _{\mathrm {L}}\gamma _{\mu }\mathrm {q}_{\mathrm {L}})], \end{aligned} \end{aligned}$$

where the subscripts \(\mathrm {L}\) and \(\mathrm {R}\) refer to the left and right chiral projections of the quark fields, respectively, and \(\eta _{\mathrm {LL}}\), \(\eta _{\mathrm {RR}}\), and \(\eta _{\mathrm {RL}}\) are taken to be 0, \({+}\,1\), or \({-}\,1\) for the different combinations that correspond to different CI models. The following CI possibilities with color-singlet couplings among quarks are investigated:

Model

\(\left( \eta _{\mathrm {LL}},\,\eta _{\mathrm {RR}},\,\eta _{\mathrm {RL}} \right) \)

\(\varLambda ^{\pm }_{\mathrm {LL}}\)

\((\pm \, 1, 0, 0)\,\)

\(\varLambda ^{\pm }_{\mathrm {RR}}\)

\((0,\pm \, 1, 0)\,\)

\(\varLambda ^{\pm }_{\mathrm {VV}}\)

\((\pm \, 1,\pm \, 1,\pm \, 1)\,\)

\(\varLambda ^{\pm }_{\mathrm {AA}}\)

\((\pm \, 1,\pm \, 1,\mp \, 1)\,\)

\(\varLambda ^{\pm }_{(\mathrm {V-A})}\)

\((0,0,\pm \, 1)\,\)

The models with positive (negative) \(\eta _{\mathrm {LL}}\) or \(\eta _{\mathrm {RR}}\) lead to destructive (constructive) interference with the QCD terms, and consequently a lower (higher) cross section, respectively. In all CI models discussed in this paper, NLO QCD corrections are employed to calculate the cross sections. In proton–proton collisions, the \(\varLambda ^{\pm }_{\mathrm {LL}}\) and \(\varLambda ^{\pm }_{\mathrm {RR}}\) models result in identical lowest order cross sections and NLO corrections, and consequently lead to the same sensitivity. For \(\varLambda ^{\pm }_{\mathrm {VV}}\) and \(\varLambda ^{\pm }_{\mathrm {AA}}\), as well as for \(\varLambda ^{\pm }_{(\mathrm {V-A})}\), the CI predictions are also identical at lowest order, but exhibit different NLO corrections and yield different sensitivities. The cijet 1.0 program [68] is used to calculate the CI terms, as well as the interference between the CI and QCD terms at NLO in QCD.

For the ADD model, two parameterizations for virtual graviton exchange are considered: Giudice–Rattazzi–Wells (GRW) [69] and Han–Lykken–Zhang (HLZ) [70]. In the GRW convention, the sum over the Kaluza–Klein graviton excitations in the effective field theory is regulated by a single cutoff parameter \(\varLambda _{\mathrm {T}}\). In the HLZ convention, the effective theory is described in terms of two parameters, the cutoff scale \(M_{\mathrm {S}}\) and the number of extra spatial dimensions \(n_{\mathrm {ED}}\). The parameters \(M_{\mathrm {S}}\) and \(n_{\mathrm {ED}}\) are directly related to \(\varLambda _{\mathrm {T}}\) [71]. We consider models with 2–6 EDs. The case of \(n_{\mathrm {ED}} =1\) is not considered since it would require an ED of the size of the radius of the solar system; the gravitational potential at such distances would be noticeably modified, and this case is therefore excluded by observation. The case of \(n_{\mathrm {ED}} =2\) is special in the sense that the relation between \(M_{\mathrm {S}}\) and \(\varLambda _{\mathrm {T}}\) also depends on the parton-parton CM energy \(\sqrt{s}\). The ADD predictions are calculated using pythia.

Quantum black hole production is studied within the framework of the ADD model, with \(n_{\mathrm {ED}} =6\) (ADD6), and the Randall–Sundrum model (RS1) [72, 73] with a single, warped extra dimension (\(n_{\mathrm {ED}} =1\)). In these models, the QBH production cross section depends on the mass of the QBH, \({M_\mathrm {Pl}}\), and the number of spatial dimensions. Since QBHs are produced with a mass threshold close to \({M_\mathrm {Pl}}\), we set the minimum QBH mass \(M_{\mathrm {QBH}}\) equal to \({M_\mathrm {Pl}}\) for simplicity. The qbh 3.0 generator [74] is used for the predictions.

To take into account the NLO QCD and EW corrections to SM dijet production when probing the ADD, QBH, and DM models, the cross section difference \(\sigma ^{\mathrm {QCD}}_{\text {NLO+EW corr}} - \sigma ^{\mathrm {QCD}}_{\text {LO}}\) is evaluated for each \(M_{\mathrm {jj}}\) and \(\chi _{\text {dijet}}\) bin and added to the SM+BSM predictions. This procedure provides an SM+BSM prediction where the QCD terms are corrected to NLO with EW corrections while the BSM terms are calculated at LO. While the ADD BSM prediction from pythia includes the interference terms of graviton exchange with QCD (obtained by subtracting the predictions \(\sigma ^{\text {ADD+QCD}}_{\text {LO}} - \sigma ^{\text {QCD}}_{\text {LO}}\)), the QBH and DM BSM predictions do not include such interference terms.

Exclusion limits on the BSM models studied in this paper are set based on the comparison of data that have not been corrected for resolution effects with both SM+BSM and SM predictions that have been folded to detector level. The comparison at detector level is done to eliminate some systematic uncertainties that are introduced during the unfolding procedure and simplifies the statistical evaluation. This procedure uses the same two-dimensional response matrix whose inverse is used for unfolding the data. It has been verified that the \(\chi _{\text {dijet}}\) distributions for SM+BSM predictions folded with the response matrix derived from SM QCD multijet predictions smeared with the double-sided Crystal Ball parameterization of the jet \(p_{\mathrm {T}}\) resolution agree with SM+BSM predictions smeared with this same parameterization. The folding procedure is equivalent to running the full detector simulation on the particle-level predictions, with the residual differences accounted for in the systematic uncertainties.

5 Systematic uncertainties

The normalized \(\chi _{\text {dijet}}\) distributions are relatively insensitive to many potential systematic effects. To present the uncertainties for the normalized shapes, the quoted values are reported for the lowest \(\chi _{\text {dijet}}\) bins, where the uncertainties and potential contributions from BSM processes are typically the largest. The main experimental uncertainty is from the jet energy scale (JES) and the main theoretical uncertainty is from the choices of \(\mu _{\mathrm {r}}\) and \(\mu _{\mathrm {f}}\) scales.

5.1 Experimental uncertainties

The overall JES uncertainty is less than 1%, and the variation of the JES as a function of pseudorapidity is less than 1% per unit \(\eta \) [47, 75] in the phase space of the analysis. The JES uncertainties related to each step in the derivation of the \(p_{\mathrm {T}}\) and \(\eta \) dependent JES corrections are taken into account independently. In this way, the correlations of the JES uncertainty sources among the \(M_{\mathrm {jj}}\) ranges and \(\chi _{\text {dijet}}\) bins are included. For the purpose of display in figures and tables, the total JES uncertainty is obtained from the quadratic sum of these uncertainty sources and is found to be 3.6% in the lowest \(M_{\mathrm {jj}}\) range and 9.2% in the highest \(M_{\mathrm {jj}}\) range.

The uncertainty from the jet \(p_{\mathrm {T}}\) resolution is evaluated by changing the width of the Gaussian core of the Crystal Ball parameterization of the response by up to ± 5% [47, 75], depending upon the jet \(\eta \), and comparing the resultant distributions before and after these changes. This uncertainty is found to be less than 1% for all \(M_{\mathrm {jj}}\). The uncertainty from the modeling of the tails of the jet \(p_{\mathrm {T}}\) resolution [76] is evaluated using a Gaussian function to parameterize the response, and we assign an uncertainty equal to half of the difference between the distributions determined from this Gaussian ansatz and the nominal correction. The size of this uncertainty is less than 1.5% for all \(M_{\mathrm {jj}}\).

Another source of uncertainty arises from the use of a parametric model to simulate the jet \(p_{\mathrm {T}}\) resolution of the detector. This uncertainty is estimated by comparing the smeared \(\chi _{\text {dijet}}\) distributions to the ones from a detailed simulation of the CMS detector using Geant4  [77], and is found to be 0.5% and 1% in the lowest and highest \(M_{\mathrm {jj}}\) ranges, respectively.

In the unfolding procedure, there is an additional systematic uncertainty introduced due to potential mismodeling of the dijet kinematic distributions in pythia. This uncertainty is evaluated using MadGraph 5_amc@nlo 2.2.2 [78] predictions, as the kinematic distributions from MadGraph 5_amc@nlo and pythia are found to bracket the data. The inverted response matrix from pythia is applied to the smeared \(\chi _{\text {dijet}}\) distributions from MadGraph 5_amc@nlo and the results are compared to the corresponding generated \(\chi _{\text {dijet}}\) distributions. The differences are observed to be less than 1.5% for all \(M_{\mathrm {jj}}\).

The effect from pileup is studied by comparing the \(\chi _{\text {dijet}}\) distributions with various numbers of pileup interactions in simulated events. The numbers are varied according to the uncertainty of the total inelastic cross section of \(\mathrm {p}\mathrm {p}\) collisions [79]. The effect on the \(\chi _{\text {dijet}}\) distributions is observed to be negligible.

5.2 Theoretical uncertainties

The uncertainties due to the choices of \(\mu _{\mathrm {f}}\) and \(\mu _{\mathrm {r}}\) scales in the NLO QCD predictions are evaluated by following the proposal in Refs. [80, 81] and changing the default choice of scales in the following 6 combinations: \((\mu _{\mathrm {f}}/\langle p_\mathrm {T} \rangle \), \(\mu _{\mathrm {r}}/\langle p_\mathrm {T} \rangle )\) = (1 / 2, 1 / 2), (1 / 2, 1), (1, 1 / 2), (2, 2), (2, 1), and (1, 2). These changes modify the predictions of the normalized \(\chi _{\text {dijet}}\) distributions by up to 8.5% and up to 19%, at small and large values of \(M_{\mathrm {jj}}\), respectively. The uncertainty in the NLO QCD predictions due to the choice of PDFs is determined from the 28 eigenvectors of CT14 using the procedure described in Ref. [82], and is found to be less than 0.2% at low \(M_{\mathrm {jj}}\) and less than 0.6% at high \(M_{\mathrm {jj}}\). The uncertainty in the strong coupling constant has a negligible impact on the normalized \(\chi _{\text {dijet}}\) distribution.

Scale and PDF uncertainties in the CI predictions are obtained using the same procedure as in the QCD predictions. In the ADD and QBH models, the scale and PDF uncertainties have a negligible impact on the limits as the signals only appear in the highest mass bins, where the statistical uncertainties dominate. The effect on the acceptance for the DM models due to the PDF uncertainty is evaluated using the 100 replica NNPDF3.0 PDF set [60] and found to be non-negligible in the \(M_{\mathrm {jj}}\) ranges with \(M_{\mathrm {jj}} >M_{\text {Med}} \) for DM mediators that have large mass and coupling. For example, for an axial-vector mediator with \(M_{\text {Med}} = 6\hbox {TeV}\) and \(g_{\mathrm {\mathrm {q}}} = 1.0\), which corresponds to a resonance with relative width of 50%, the uncertainty is 14% in the \(M_{\mathrm {jj}} >\!6.0\hbox {TeV}\) bin.

Although the uncertainties are treated separately in the statistical analysis of the data, for display purposes in tables and figures we calculate the total experimental and theoretical uncertainty as the quadratic sum of the contributions due to the JES, the jet \(p_{\mathrm {T}}\) resolution, the modeling of both the detector response and the dijet kinematics, and the contributions from \(\mu _{\mathrm {f}}\), \(\mu _{\mathrm {r}}\), and the PDFs. A summary of the leading experimental systematic uncertainties is provided in Table 1. The theoretical uncertainties quoted in the table apply to the QCD prediction. As shown in Table 1, systematic uncertainties dominate the total uncertainty in low \(M_{\mathrm {jj}}\) regions, while the statistical uncertainty dominates in high \(M_{\mathrm {jj}}\) regions.

Table 1 Summary of the leading experimental and theoretical uncertainties in the normalized \(\chi _{\text {dijet}}\) distributions, in percent. While the statistical analysis represents each uncertainty through a change in the \(\chi _{\text {dijet}}\) distribution correlated among all \(\chi _{\text {dijet}}\) bins, this table summarizes each uncertainty by a representative value to show their relative contributions. For the lowest and highest dijet mass bins, the relative shift is given for the lowest \(\chi _{\text {dijet}}\) bin. In the highest dijet mass bin, the dominant experimental contribution corresponds to the statistical uncertainty, while the dominant theoretical contribution is from the uncertainty in scales

6 Results

In Figs. 1 and 2 the measured normalized \(\chi _{\text {dijet}}\) distributions for all mass bins unfolded to particle level are compared to NLO predictions with EW corrections. No significant deviation from the SM prediction is observed. The distributions are also compared to predictions for QCD+CI with CI scales equal to 14TeV, QCD+ADD with \(\varLambda _{\mathrm {T}}\ (\mathrm {GRW})\ =10\hbox {TeV}\), QCD+QBH with \(M_{\mathrm {QBH}} ~(\mathrm {ADD6}) = 8\hbox {TeV}\), and QCD+DM with \(M_{\mathrm {Med}}=2\), 3 and 5TeVand \(g_{\mathrm {\mathrm {q}}} =1.0\). The signal distributions are only shown for the \(M_{\mathrm {jj}}\) ranges that contribute to the sensitivity for the BSM searches.

Fig. 1
figure 1

Normalized \(\chi _{\text {dijet}}\) distributions in the three highest mass bins. Unfolded data are compared to NLO predictions (black dotted line). The error bars represent statistical and experimental systematic uncertainties combined in quadrature. The ticks on the error bars correspond to the experimental systematic uncertainties only. Theoretical uncertainties are indicated as a gray band. Also shown are the predictions for various CI, ADD, QBH, and DM scenarios. The lower panels show the ratio of the unfolded data distributions and NLO predictions

Fig. 2
figure 2

Normalized \(\chi _{\text {dijet}}\) distributions in the four lower mass bins. Unfolded data are compared to NLO predictions (black dotted line). The error bars represent statistical and experimental systematic uncertainties combined in quadrature. The ticks on the error bars correspond to the experimental systematic uncertainties only. Theoretical uncertainties are indicated as a gray band. Also shown are the predictions for various CI, ADD, and DM scenarios. The lower panels show the ratio of the unfolded data distributions and NLO predictions

The asymptotic approximation [83] of the \(\mathrm {\text {CL} _s}\) criterion [84, 85] is used to set exclusion limits on the parameters for the BSM models [86]. The limits obtained using this approximation were tested against the \(\mathrm {\text {CL} _s}\) limits obtained using ensembles of pseudo experiments for several of the models examined, and the differences were found to be negligible. The likelihoods \(L_\mathrm {QCD}\) and \(L_\mathrm {QCD+BSM}\) are defined for the respective QCD-only and QCD+BSM hypotheses as a product of Poisson likelihood functions for each bin in \(\chi _{\text {dijet}}\). The predictions for each \(M_{\mathrm {jj}}\) range are normalized to the number of observed events in that range. Systematic uncertainties are treated as nuisance parameters in the likelihood model. Following Ref. [17], the nuisance parameters are profiled with respect to the QCD-only and QCD+BSM models by maximizing the corresponding likelihood functions. The p-values for the two hypotheses, \(P_\mathrm {QCD+BSM}(q \ge q_\text {obs})\) and \(P_\mathrm {QCD}(q \le q_\text {obs})\), are evaluated for the profile log-likelihood ratio \(q = -\,2 \ln ({L_\text {QCD+BSM}}/{L_\mathrm {QCD}})\). Limits on the QCD+BSM models are set based on the quantity \(\mathrm {\text {CL}}_\mathrm {s} = P_\mathrm {QCD+BSM}(q\ge q_\text {obs}) / (1-P_\mathrm {QCD}(q \le q_\text {obs}))\), which is required to be less than 0.05 for a 95% confidence level (\(\text {CL}\)) of exclusion. Because of the large number of events in the low-\(M_{\mathrm {jj}}\) range, which constrain the systematic uncertainties, we obtain 2–30% better observed limits on the BSM scales and masses compared to the limits obtained using the method in the predecessor of this search reported in Ref. [29], where the nuisance parameters were marginalized rather than profiled.

In the limit calculations, not all \(M_{\mathrm {jj}}\) ranges are included in the likelihoods; only those that improve the expected limits by more than 1% are used. We use mass bins with \(M_{\mathrm {jj}} >3.6\hbox {TeV}\) for the CI models, \(M_{\mathrm {jj}} >4.2\hbox {TeV}\) for the ADD models, and \(M_{\mathrm {jj}} >4.8\hbox {TeV}\) for the QBH models. For the DM mediators, we use mass bins that cover the \(M_{\mathrm {jj}}\) range of 0.5\(M_{\text {Med}}\)–1.2\(M_{\text {Med}}\). The exclusion limits on the BSM models are determined using detector-level \(\chi _{\text {dijet}}\) distributions and theoretical predictions at detector level. By using the detector-level \(\chi _{\text {dijet}}\) distributions, each bin of the \(\chi _{\text {dijet}}\) distributions can be modeled by a Poisson likelihood function, while at particle level, the unfolded data distributions have correlations among the dijet mass bins. As a cross-check, the limits are also determined for the case where the unfolded \(\chi _{\text {dijet}}\) distributions, approximated by Poisson likelihood functions, and particle-level theoretical predictions are used in the limit extraction procedure. The resulting observed limits on the BSM scales and masses are found to be more stringent than those determined at detector level by 1–10%, depending on the model. The agreement of the data with QCD predictions is quantified by calculating \(P_\text {QCD}(q < q_\text {obs})\) for each mass bin separately. The largest excess is found in the first data point of the > 6.0TeVmass bin, with a significance of 1.8 standard deviations. When combining mass bins in the various QCD+BSM models under study, the largest significances are found to be 2.7–2.8 standard deviations for the QCD+DM model with \(M_{\text {Med}} = 4.5\)–6TeVand \(g_{\mathrm {\mathrm {q}}} =1.0\).

Figure 3 shows the 95% \(\text {CL}\) upper limits on \(g_{\mathrm {\mathrm {q}}}\) as a function of the mass of the vector or axial-vector DM mediator with \(g_{\mathrm {DM}} =1.0\) and \(m_{\mathrm {DM}} =1\hbox {GeV}\). The corresponding limits on the width of the mediators are shown on the vertical axis on the right-hand side of Fig. 3. The degradation of the limits below \(M_{\text {Med}} =2.5\hbox {TeV}\) and above \(M_{\text {Med}} =4\hbox {TeV}\) can be explained as follows. For resonance masses below the lower \(M_{\mathrm {jj}}\) boundary of the analysis at 2.4TeV, the acceptance increases rapidly as a function of resonance mass (e.g., from 1.4% at \(M_{\text {Med}} =2\hbox {TeV}\) to 16% at \(M_{\text {Med}} =2.5\hbox {TeV}\), for \(g_{\mathrm {\mathrm {q}}} =0.5\)), resulting in the improvement of the limit on \(g_{\mathrm {\mathrm {q}}}\) as a function of resonance mass. For large values of resonance mass and width (e.g., for \(M_{\text {Med}} >4\hbox {TeV}\) and \(g_{\mathrm {\mathrm {q}}} >0.5\)), the mediator is primarily produced off-shell with a mass less than the \(M_{\mathrm {jj}}\) boundary of the analysis at 2.4TeV. The acceptance for high resonance masses thus decreases as a function of resonance width (e.g., for \(M_{\text {Med}} =5\hbox {TeV}\), from 25% at \(g_{\mathrm {\mathrm {q}}} =0.5\) to 8% at \(g_{\mathrm {\mathrm {q}}} =1.5\)), resulting in the fast deterioration of the limit on \(g_{\mathrm {\mathrm {q}}}\) at high resonance masses. The observed limit above 5TeVis at \(\varGamma /M_{\mathrm {Med}} \ge 1\), thus in a region where the simplified model of a mediator particle is no longer valid. For \(M_{\text {Med}}\) between 2.0 and 4.6TeV, this search excludes couplings \(1.0 \le g_{\mathrm {\mathrm {q}}} \le 1.4\), which are not accessible via dijet resonance searches.

The limits for \(M_{\text {Med}}\) at arbitrary \(m_{\mathrm {DM}}\) and \(g_{\mathrm {DM}}\) can be calculated based on the fact that at fixed mediator production cross sections, changes in the width of the DM decay channel will lead to changes in the width of the quark decay channel. For the models with \(g_{\mathrm {\mathrm {q}}} =1.0\), \(g_{\mathrm {DM}} =1.0\), and \(2m_{\mathrm {DM}} <M_{\text {Med}} \), in which the total width of the mediator is dominated by the width of the quark decay channel due to the large number of possible quark flavors and colors, the exclusion range for \(M_{\text {Med}}\) has little dependence on \(m_{\mathrm {DM}}\). For the models with \(2m_{\mathrm {DM}} \ge M_{\text {Med}} \), the width of the DM decay channel is assumed to be zero. The resulting exclusion regions for vector and axial-vector mediators with \(g_{\mathrm {\mathrm {q}}} =1.0\) and \(g_{\mathrm {DM}} =1.0\) in the \(m_{\mathrm {DM}}\) and \(M_{\text {Med}}\) plane are shown in Fig. 4.

Fig. 3
figure 3

The 95% \(\text {CL}\) upper limits on the quark coupling \(g_{\mathrm {\mathrm {q}}}\), as a function of mass, for an axial-vector or vector DM mediator with \(g_{\mathrm {DM}} =1.0\) and \(m_{\mathrm {DM}} =1\hbox {GeV}\). The observed limits (solid), expected limits (dashed) and the variation in the expected limit at the 1 and 2 standard deviation levels (shaded bands) are shown. A dotted horizontal line shows the coupling strength for a benchmark DM mediator with \(g_{\mathrm {\mathrm {q}}} =1.0\). The corresponding limits on the width of the mediators are shown on the vertical axis on the right-hand side of the figure

Fig. 4
figure 4

The 95% \(\text {CL}\) observed (red) and expected (blue) excluded regions in the plane of \(m_{\mathrm {DM}}\) and \(M_{\text {Med}}\), for a vector mediator (upper) and an axial-vector mediator (lower) for a DM benchmark model with \(g_{\mathrm {DM}} =g_{\mathrm {\mathrm {q}}} =1.0\). These are compared to constraints from the cosmological relic density of DM (gray) determined from astrophysical measurements [87], using MadDM. In the hatched area, DM is overabundant. The observed and expected lower bounds for \(M_{\text {Med}}\) overlap with each other

The observed and expected exclusion limits at 95% \(\text {CL}\) on different CI, ED, QBH, and DM models obtained in this analysis are listed in Table 2. The observed limits are less stringent than the expected limits because of the upward fluctuation in the measured distributions compared to the theoretical predictions. The limits on all models are more stringent than those obtained from data collected by CMS in 2015 [29].

Table 2 Observed and expected exclusion limits at 95% \(\text {CL}\) for various CI, ADD, QBH, and DM models. The 68% ranges of expectation for the expected limit are given as well. For the DM vector mediator, couplings \(g_{\mathrm {DM}} =1.0\), \(g_{\mathrm {\mathrm {q}}} \ge 1\) and a DM mass of 1GeVare assumed and a range of masses instead of a lower limit is quoted

7 Summary

A search has been presented for physics beyond the standard model, based on normalized dijet angular distributions obtained in 2016 from proton–proton collisions at the LHC. The data sample corresponds to an integrated luminosity of 35.9\(\,\text {fb}^{-1}\). The angular distributions, measured over a wide range of dijet invariant masses, are found to be in agreement with the predictions of perturbative quantum chromodynamics. The results are used to set 95% confidence level lower limits on the contact interaction scale for a variety of quark compositeness models, the ultraviolet cutoff in models of extra spatial dimensions, the minimum mass of quantum black holes, and the mass and couplings in dark matter models. For the first time, lower limits between 2.0 and 4.6TeVare set on the mass of a dark matter mediator for (axial-)vector mediators, for the universal quark coupling \(1.0\le g_{\mathrm {\mathrm {q}}} \le 1.4\). This region is not accessible through dijet resonance searches. The lower limits for the contact interaction scale \(\varLambda \) range from 9.2 to 22.4TeV. The lower limits on the ultraviolet cutoff in the Arkani–Hamed–Dimopoulos–Dvali model are in the range of 8.5–12TeV, and are the most stringent limits available. Quantum black hole masses below 8.2TeVare excluded in the model with six large extra spatial dimensions, and below 5.9TeVin the Randall–Sundrum model with a single, warped extra dimension. To facilitate comparisons with the predictions of other models, the angular distributions, corrected to particle level, are published in HEPData.