Search for resonant and nonresonant production of pairs of dijet resonances in proton-proton collisions at $\sqrt{s}$ = 13 TeV

A search for pairs of dijet resonances with the same mass is conducted in final states with at least four jets. Results are presented separately for the case where the four jet production proceeds via an intermediate resonant state and for nonresonant production. The search uses a data sample corresponding to an integrated luminosity of 138 fb$^{-1}$ collected by the CMS detector in proton-proton collisions at $\sqrt{s}$ = 13 TeV. Model-independent limits, at 95% confidence level, are reported on the production cross section of four-jet and dijet resonances. These first LHC limits on resonant pair production of dijet resonances via high mass intermediate states are applied to a signal model of diquarks that decay into pairs of vector-like quarks, excluding diquark masses below 7.6 TeV for a particular model scenario. There are two events in the tails of the distributions, each with a four-jet mass of 8 TeV and an average dijet mass of 2 TeV, resulting in local and global significances of 3.9 and 1.6 standard deviations, respectively, if interpreted as a signal. The nonresonant search excludes pair production of top squarks with masses between 0.50 TeV to 0.77 TeV, with the exception of a small interval between 0.52 and 0.58 TeV, for supersymmetric $R$-parity-violating decays to quark pairs, significantly extending previous limits. Here, the most significant excess above the predicted background occurs at an average dijet mass of 0.95 TeV, for which the local and global significances are 3.6 and 2.5 standard deviations, respectively.


Introduction
The potential to discover physics beyond the standard model (BSM) is the motivation of searches for dijet resonances. These searches look for a pair of jets coming from the decay of a new particle, X.
The simplest searches, which have been conducted at the CERN LHC [1-18] using strategies reviewed in Ref. [19], are for s-channel single production of dijet resonances produced in proton-proton (pp) collisions, where there are two jets in the final state: pp → X → jj. Because of this simple topology, where X is produced directly from the partons in the proton and must decay to the same partons, searches for single production place tight constraints on many models of BSM physics.
A more complex process is pair production of a dijet resonance X decaying to two jets and therefore yielding four jets in the final state. Two modes of production of the pair of new particles are possible, resonant production with an intermediate resonance Y and non-resonant production, both shown in Fig. 1, and each mode presents a unique opportunity for discovery. We present a generic search for pairs of narrow dijet resonances in each production mode with data corresponding to an integrated luminosity of 138 fb −1 collected in 2016-2018 with the CMS detector at the LHC, and we use results obtained from each mode to constrain a benchmark model, as discussed below. Prior searches were for Lorentz-boosted low mass resonance pairs, produced via both resonant [20] and nonresonant [21] mechanisms, and have been conducted using jet substructure techniques to identify boosted resonances inside of jets. This analysis considers pairs of resolved dijet resonances, X, where both jets within each dijet resonance are individually reconstructed, allowing the resonant and nonresonant searches to be sensitive to high resonance masses: X masses greater than 0.5 TeV and Y masses greater than 2.0 TeV. First, we consider resonant production of pairs of dijet resonances, where the intermediate state is a massive new particle, Y, decaying to identical dijet resonances, X. As a benchmark, we consider a diquark model [22], where the intermediate state is a diquark, S, produced from the annihilation of two up quarks, and the dijet resonance is a vector-like quark, χ, that decays to an up quark (u) and a gluon (g) only; uu → S → χχ → (ug)(ug).
This scalar diquark is a good benchmark because it is produced with a large cross section, due to the high probability of finding up quarks at high fractional momentum within the proton. This is the first time CMS has searched for resonant pair production of dijet resonances via high mass intermediate states. We are further motivated to conduct this search by an event, a clear candidate for resonant production of pairs of dijet resonances, that was seen and presented by a CMS search in the dijet final state [1]. As discussed in Ref. [22], that search was not optimized for this signal and could not quantify its significance, motivating a dedicated search in the four-jet final state. Here we present an optimized search for resonant production of pairs of dijet resonances, exploring a wide range of four-jet and dijet masses.
As a benchmark, we consider the R-parity-violating (RPV) supersymmetry (SUSY) model [23], where the dijet resonance is a top squark, t, with RPV decays to down (d) and strange (s) quarks pp → t t * → (ds)(ds). Tabulated results are provided in the HEPData record for this analysis [33].

The CMS detector
A detailed description of the CMS detector and its coordinate system, including definitions of the azimuthal angle φ and pseudorapidity η, is given in Ref. [34]. The central feature of the CMS apparatus is a superconducting solenoid of 6 m internal diameter providing an axial magnetic field of 3.8 T. Within the solenoid volume are located the silicon pixel and strip tracker, and the barrel and endcap calorimeters (|η| < 3.0), where these latter detectors consist of a lead tungstate crystal electromagnetic calorimeter, and a brass and scintillator hadron calorimeter. An iron and quartz-fiber hadron calorimeter is located in the forward region (3.0 < |η| < 5.0), outside the solenoid volume. The muon detection system covers |η| < 2.4 with up to four layers of gas-ionization detectors installed outside the solenoid and embedded in the layers of the steel flux-return yoke.

Simulated data samples
A background sample of quantum chromodynamics (QCD) multijet events is produced, for the optimization of the search and qualitative comparisons with the observed data. The generation starts from the leading order (LO) QCD 2 → 2 processes of jet production, and includes additional jets from QCD initial-and final-state radiation within the parton shower. These background predictions are produced with the PYTHIA 8.205 program, with the CUETP8M1 tune [35,36], using the parton distribution function (PDF) set NNPDF2.3LO [37].
The benchmark model for the resonant search is based on the aforementioned model of a diquark decaying to pairs of vector-like quarks, defined in Eq. (2). Specifically, we consider a scalar diquark resonance S which decays to two vector-like quarks χ, which each then decays to a quark and gluon pair: S → χχ → (ug)(ug). The MADGRAPH5 aMC@NLO 2.6.5 [38] generator, with additional code specifying the model [22], is used to generate these events and MADSPIN is used for the vector-like quark decay to a quark and a gluon. Events for diquark masses between 2 and 9 TeV and for vector-like quark masses between 0.2 and 3.8 TeV are generated.
The benchmark model for the nonresonant search is the pair production of top squarks defined in Eq. (3). In the RPV SUSY model [39], the top squarks decay to d and s quarks via a threegeneration coupling λ 312 with 100% branching fraction: B(ds) = 1. For this model, we also use MADGRAPH5 aMC@NLO to generate events with up to one additional jet using the MLM matching procedure [40]. Samples with top squark masses between 0.5 and 3.0 TeV are used.
For both benchmark signal simulations, PYTHIA 8.205 [41] is used to produce the parton shower and the resulting final-state particles. For the diquark simulated samples the CP2 underlying event tune [35] is used with the NNPDF3.1LO PDF set [42], and for the top squark samples the CP5 tune [43] is used with the NNPDF3.1NNLO PDF set [42], a next-to next-to leading order PDF.
The simulation of the CMS detector for all samples is handled by GEANT4 [44]. All samples include the effects of pileup, additional pp interactions in the same or adjacent bunch crossings. This pileup distribution in simulation is weighted to match the one observed in data.

Event reconstruction and selection
A particle-flow (PF) event algorithm aims to reconstruct and identify each individual particle in an event, with an optimized combination of information from the various elements of the CMS detector [45]. Particles are classified as muons, electrons, photons, charged or neutral hadrons. To reconstruct jets from the input particles, the anti-k T algorithm [46,47] is used with a distance parameter of 0.4 as implemented in the FASTJET package [48]. The primary vertex (PV) is taken to be the vertex corresponding to the hardest scattering in the event, evaluated using tracking information alone, as described in Section 9.4.1 of Ref. [49]. Charged PF candidates not originating from the primary vertex are removed prior to the jet finding.
Events are selected using a two-tier trigger system [50]. Events satisfying loose jet requirements at the first-level (L1) trigger are examined by the high-level trigger (HLT) system. At L1, singlejet triggers that require at least one jet in the event to exceed a predefined transverse momentum (p T ) threshold are used. Triggers that require H T to exceed a threshold, where H T is the scalar sum of jet p T for all jets in the event with p T > 30 GeV and |η| < 3.0, are also used. The HLT requires H T > 1050 GeV or at least one jet reconstructed with an anti-k T distance parameter of 0.8 and p T > 550 GeV.
The jet momenta and energies are corrected using calibration factors obtained from simulation, test beam results, and pp collision data at √ s = 13 TeV. The methods described in Ref. [51] are used and all in-situ calibrations are obtained from the current data. Jets are required to have p T > 80 GeV and |η| < 2.5. The four jets with the largest p T are defined as the four leading jets. Jet identification criteria are applied to remove spurious jets associated with calorimeter noise, as well as those associated with muon and electron candidates that are either misreconstructed or isolated [52]. For all jets, we require that the component of the jet energy arising from neutral hadrons and the component arising from photons, are each less than 90% of the total jet energy. For jets within the fiducial tracker coverage we additionally require the jet to have nonzero charged-hadron energy, and electron and muon energies to be less than 90% and 80% of the total jet energy respectively. An event is rejected if any of the four leading jets fails these jet identification criteria.
The dijet pairs are constructed from the four leading jets. In the following definitions the sub-scripts 1, and 2 refer to the two dijet systems: dijet 1 and dijet 2. There are three possible combinations of four jets making two dijet systems. A jet pairing algorithm is used to choose one of these three combinations to form the pair of dijets, and the jets are relabeled with indices a, b, c and d. The pairing algorithm uses the η-φ space separations between the jets in each dijet, ∆R 1 and ∆R 2 , given by where indices j a and j b refer to jets in the first dijet system (subscript 1), and indices j c and j d refer to jets in the second dijet system (subscript 2). We find that the optimal pairing algorithm is the one minimizing the following measure of the η-φ space separations of the jets in the event: This pairing algorithm was first used in the search for nonresonant production of pairs of dijet resonances [21]. It is motivated by the expectation that the jets from a dijet resonance will be closer together than uncorrelated jets, minimizing the separation of the jets, while the offsets of 0.8 in Eq. (6) reduce pairings where the jets overlap in η-φ space.
Once the jet pairing is complete, events are selected that satisfy: These requirements reject background from hard multijet processes produced by QCD, which does not naturally give small separations in η-φ space between the jets in each dijet in a four-jet event. In addition, we require the η separation of the two dijets to satisfy: which further suppresses background from QCD t-channel production compared to the signal from s-channel production. Finally, the asymmetry in dijet mass between the two dijets is required to be small, selecting dijet pairs of approximately equal mass, a property of pairs of identical resonances, but not of the QCD background. Here, m 1 and m 2 denote the reconstructed mass of each dijet pair. Throughout the paper, the small letter m denotes reconstructed dijet and four-jet masses, and the capital letter M denotes true resonance masses.
An analysis optimization procedure was performed using Monte Carlo (MC) simulations of signals and background. In addition to the optimal pairing method that minimized ∆R in Eq. (6), we also considered a second method that minimized the asymmetry, a third utilizing both ∆R and asymmetry, a fourth minimizing the difference between ∆R 1 and ∆R 2 , a fifth minimizing the sum of the average dijet mass and the difference in dijet mass, and finally a sixth method minimizing both the average dijet mass and the difference in dijet mass. We varied the thresholds for the event selection requirements in Eqs. (5-9), and determined the expected signal significance at both low and high resonance masses, for both resonant and nonresonant signals. The optimal dijet pairing method and event selections, presented in Eqs. (6-9), are the ones that maximized the expected signal significance for the full analysis. The different jet pairings explored do not significantly affect the resonant search, whose final observable is the four-jet mass, and give similar results for the non-resonant search.

Analysis strategy
Natural variables for the two searches are derived from the properties of the production processes. The two dijet resonances XX have equal mass, so one-dimensional (1D) distributions of the average dijet mass m jj = (m 1 + m 2 )/2 can be used to search for the nonresonant production of this pair. Resonant pair production, Y → XX, additionally produces a localized excess in a second variable, the four-jet mass m 4j , clearly visible in the two-dimensional (2D) distribution of m jj vs. m 4j . Figure 2 shows the 2D distribution of the observed data for m 4j vs. m jj . Near m 4j ∼ 8 TeV and m jj ∼ 2 TeV we observe two events, candidates for a resonant signal. For comparison, Fig. 2 also shows the 68 and 95% probability contours of the signal, respectively, from a simulation of our benchmark model of resonant pair production, for a diquark of mass 8.4 TeV and a vectorlike quark mass of 2.1 TeV. The challenge in evaluating the significance of this signal, or any signal, is to understand the background. Perturbative QCD predictions of multijet production have too many theoretical uncertainties to reliably model this background. We instead use an approach that does not rely on simulation, by fitting a smooth background parameterization function to the observed data. We do this with a set of 1D distributions that span the 2D space.  The seemingly natural approach of binning the search would be to construct the set of 1D distributions of the observed number of events as a function of m 4j , in bins of m jj , or vice versa. However, Fig. 2 shows that due to correlations between those two variables this approach would introduce sculpting of the distributions, which would limit the range of new particle masses we could search for. Hence, we define a new dimensionless variable, α = m jj /m 4j , which is the ratio of the two defined variables of resonant production, and a measure of the boost of the two dijet systems for both production modes. For fixed values of α, the number of events observed decreases smoothly with increasing m 4j , or increasing m jj . Therefore, binning in a α provides a wide unbiased range of m 4j or m jj because the distributions are not sculpted or truncated by the boundaries. Thus α is a good variable to bin the 2D data to construct a set of unbiased 1D distributions for the data-driven search. In addition, although the dijet pairs from a nonresonant signal are not produced at a fixed value of α, using a binning in the nonresonant search exploits the signal distribution in α. As can be seen in Fig. 4, properly reconstructed signals (i.e., pairs of dijets with average masses near the generated t mass) tend to be produced with high α. This is especially true for higher-mass t pairs, and it is because the lower values of α only occur when the individual ts are produced with high p T . Hence, examining the four-jet and average dijet mass distributions in bins of α is expected to significantly increase the sensitivity of both searches.
The average dijet mass and the four-jet mass are the variables used for the signal extraction in the nonresonant and resonant search respectively, as described in more detail below. We then optimize, for maximal signal significance, both the number of α bins considered and the final observable used to search for localized excesses. The resonant search consists of thirteen simultaneous 1D searches using the m 4j distributions within the thirteen α bins shown in Fig. 3, and similarly the nonresonant search uses the m jj distributions in the three α bins in Fig. 4. For the resonant search the lowest α value we consider is 0.1 with a bin width of 0.02, which is comparable to the reconstructed α resolution, up to an α of 0.34, and then we utilize an overflow bin for all higher α values. For the nonresonant search, the optimization yielded a bin width of 0.1, for two bins between α of 0.15 and 0.35, and again a final bin is used for higher α values.
Lastly, the searches are performed above mass thresholds for which the trigger is found to be fully efficient. In the case of the resonant analysis, the search is performed for four-jet masses m 4j > 1.6 TeV for all α values. In the case of the nonresonant analysis, the thresholds in average dijet mass vary with α: m jj is required to be greater than 0.35, 0.53, and 0.61 TeV for the three bins with lower α edges of 0.15, 0.25, and 0.35, respectively. In Fig. 4, this requirement on m jj is removed for illustrative purposes, and the same events as in Fig. 3 are shown as a function of α and m jj . Figure 5 shows the shapes of signals in the resonant and nonresonant searches for various S and top squark masses, respectively, for the α bin with the largest signal present. The distributions of four-jet mass, for signals in the resonant search, peak just below the input generator mass and have a long tail to low masses from radiation. The distributions of average dijet mass, for dijet resonances in the nonresonant search, also have a peak just below the input generator mass from the correct pairings of the four jets, and an additional peak at lower masses from the incorrect pairings.

Signal simulations
We estimate the pairing efficiency by measuring the fraction of events in which the average reconstructed dijet mass agrees with the true dijet resonance mass within a window of ±3 sigma, where sigma is the average experimental resolution of the reconstructed dijet mass. Normalized yield/TeV  Efficiency × Acceptance λ RPV coupling Figure 6: The products of acceptance and efficiency of a resonant signal vs. the diquark mass (left), and a nonresonant signal vs. the top squark mass (right), inclusively, i.e., for all α values, and for the three α bins that contain the majority (≥85%) of the signal. The case when the efficiency of the mass selection is unity is also shown as a solid line.
as the fraction of events passing the kinematic selection criteria. The efficiency is the fraction of signal events satisfying the mass threshold previously discussed: four-jet mass requirement for the resonant search, dijet mass selection for the nonresonant search. Figure 6 also shows the signal acceptance alone, the curve where the efficiency is equal to 1. Cross section upper limits in this paper are quoted within this acceptance and corrected for efficiency. The highest signal acceptance for the resonant search, with α true = M(X)/M(Y) = 0.25, is naturally for the corresponding bin with 0.24 < α < 0.26. For the nonresonant search, which does not have an α true or a natural value of α, the highest signal acceptance is for the bin with 0.25 < α < 0.35. The acceptance, which originates from the dijet matching and angular separation requirements of Eqs. (5-9), is approximately independent of the resonance mass in both searches. However, there is a small decrease in the acceptance at the very highest values of signal mass in Fig. 6 for both searches, which is caused by the mass asymmetry requirement. That requirement is more likely to reject events with incorrect dijet pairings, which are more frequent at high signal masses, as previously discussed.

Background estimation
The background in this search is derived exclusively from data. As in previous resonance searches [3, 5-18, 53], we fit empirical functional forms to the observed multijet mass distributions. Our primary background function, which gives the best fit to the differential cross section data without including any signal component, is the modified dijet function (ModDijet-3p) given by dσ/dm = p 0 (1 − x 1/3 ) p 1 /x p 2 , where p 0 , p 1 , p 2 are free parameters. Here, x = m/ √ s, with m being the mass used for the search, which is the four-jet mass in the resonant search and the average dijet mass in the nonresonant search. We use two alternate background functions to estimate the systematic uncertainty in the background from variations in the functional form: the dijet function (Dijet-3p) given by dσ/dm = p 0 (1 − x) p 1 /x p 2 , and the power-law times exponential function (PowExp-3p) given by dσ/dm = p 0 e −p 1 x /x p 2 . Figure 7 shows the three functional forms, fitted to examples of the observed differential cross section, as functions of four-jet mass, and the simulation of the QCD background, for three out of the thirteen α bins in the resonant search. All functions give a good fit to the data, and a similar prediction of the background to the LO QCD simulation, and none of these background estimates appear to account for the two events observed at a four-jet mass of 8 TeV. The three α bins shown in Fig. 7 contain 85% of the signal for a four-jet resonance, Y, decaying to a pair of dijet resonances, X, with α true = M(X)/M(Y) = 0.25, for which representative signal shapes are also shown in Fig. 7. Figure 8 shows the pulls from the fit with the ModDijet-3p function for all thirteen α bins, along with the χ 2 and number of degrees of freedom (NDF) of each individual fit. All the fits adequately describe the data, with individual p-value [54] ranging from 0.14 to 0.91. The combined p-value is 0.18 for a simultaneous fit of the data with the background function, including all thirteen α bins. This came from conducting a goodness-offit test on pseudo-experiments, utilizing a binned likelihood ratio with respect to the saturated model as the test statistic. This probability, and the pulls of the individual fits, indicate that the background model gives a good description of the data overall.
For the nonresonant search, Fig. 9 shows the data, fits, QCD background simulation, and signal shapes for all three α bins. The p-values of the background fits for the individual α bins range from 0.02 to 0.87, and the combined p-value, calculated with the same statistical procedure as in the resonant search, is 0.04. Here the low p-values originate from the large absolute values of the pulls, near an average dijet mass of roughly 1 TeV for 0.25 < α < 0.35 in Fig. 9. Similar pulls can also be seen in the corresponding α bins of Fig. 8, near a four-jet mass of roughly 3-4 TeV. Only the lowest α bin in the nonresonant search extends to sufficiently small values of average dijet mass to contain a noticeable peak for a 0.6 TeV top squark, as shown in Fig. 9. It is the only α bin sensitive to top squark masses between 0.5 and 0.6 TeV.   Figure 9: The average dijet mass distributions for the data (points), within the three α bins of the nonresonant search, compared with LO QCD (green histogram) and fitted with three functions: a power-law times an exponential (red dotted), the dijet function (red dashed), and the modified dijet function (red solid), each function with three free parameters. Examples of predicted top squark pair production signals are shown, with cross sections equal to the expected SUSY cross sections, for top squark masses of 0.6 (blue solid), 1.0 (blue dotted), and 2.0 TeV (blue dashed). The lower panels show the pulls from the fit of the modified dijet function to the dijet mass distribution, calculated using the statistical uncertainty of the data, the aforementioned signals, and the goodness of fit measure χ 2 /NDF, for all three α bins of the nonresonant search.

Results
The search uses a fit of the background function plus the simulated signal shape to the data, taking into account statistical and systematic uncertainties. As discussed, the data are the fourjet mass distributions in thirteen bins of α for the resonant search, and the average dijet mass distributions in three bins of α for the nonresonant search. The fit is done simultaneously for all the bins of α within each search. The results of the search are the expected and observed upper limits on the signal cross section for all experimentally accessible values of resonance mass. For some values of resonance mass, the best fit cross section and observed statistical significance of the signal hypothesis are also of interest.
The dominant sources of systematic uncertainty in the signal modeling are the jet energy scale, the jet energy resolution, and the integrated luminosity. The uncertainties with the largest impact on the extracted signal yield are the ones relating to background, with the signal modeling related ones being negligible. The uncertainty in the jet energy scale is below 2% for all values of the four-jet and dijet masses. This uncertainty is propagated to the limits by shifting the four-jet and dijet mass shapes for the signals by ±2%. The uncertainty in the jet energy resolution translates into an uncertainty of 10% in the resolution of the four-jet and dijet masses [51], and is propagated to the limits by observing the effect of increasing the reconstructed width of the four-jet and dijet mass shapes for the signals by 10%. The total integrated luminosity has an uncertainty of 1.6%, the improvement in precision relative to Refs. [55][56][57] reflects the uncorrelated time evolution of some systematic effects, and is propagated to the normalization of the signal. Changes in the background function used, and the values of the background function parameters as estimated from the fit, can introduce changes in the signal yield. All these systematic uncertainties are included in the extracted cross section, limits and significance, as discussed in the next paragraph.
For the signal extraction, we use a multibin counting-experiment likelihood, which is a product of Poisson distributions corresponding to different bins. We evaluate the likelihood independently at each value of resonance mass from 1.6 to 10 TeV in 0.1 TeV steps for the resonant search and, for the nonresonant search, in 0.025 TeV steps between 0.5 and 1.0 TeV, 0.05 TeV steps between 1.0 and 1.8 TeV, and 0.2 TeV steps between 1.8 and 3.0 TeV. These step sizes are comparable to the experimental mass resolution at the lower edge of each resonance mass interval. The sources of systematic uncertainty are implemented as nuisance parameters in the likelihood model, with Gaussian constraints for the jet energy scale and resolution, and lognormal constraints for the integrated luminosity. The main source of the uncertainty in the background modeling comes from the choice of the background function, and the range of possible values of the parameters of the background function. The parameters of the empirical functional form used to describe the background are considered as freely floating nuisances, and are evaluated via profiling. The discrete profiling method [58] is used for considering the choice of the functional form as a discrete nuisance parameter, which is profiled in an analogous way to continuous nuisance parameters.
Next, we proceed in setting limits, and estimating significances. The modified frequentist criterion [59,60] is used to set upper limits on the signal cross sections, following the prescription described in Refs. [61,62], using the asymptotic approximation of the test statistic. In cases where the small number of events would make this approximation less accurate, the limit and significance estimation with respect to the background-only hypothesis is performed using test-statistic distributions derived from pseudo-experiments. Figures 10-12 show the model-independent observed and expected upper limits at 95% confidence level (CL) on σ B A, which is the product of the cross section (σ), the branching fraction (B), and the acceptance (A) for the kinematic requirements. Figure 10 shows the limits on resonant production for twelve values of α true , chosen at the center of each bin of α, and Fig. 11 shows the limit for α true = 0.25, for which the maximum significance is observed. All upper limits presented can be compared to the predictions of σ B A to determine mass limits on new particles. In Figs. 10 and 11 we compare these limits from the resonant search, to the predicted cross section times branching fraction for the scalar diquark model, multiplied by the acceptance in Fig. 6. For the prediction we use a next-to-leading order (NLO) calculation with NNPDF NLO 3.1 PDFs, discussed in Ref. [22], for diquark masses greater than 6 TeV. For lower diquark masses the NLO calculation is not available, so we use a MADGRAPH5 aMC@NLO calculation multiplied by a K factor of 1.47, which is the ratio of the NLO and MADGRAPH cross sections at 6 TeV. In Fig. 12, the top squark signal cross section is calculated from a simplified topology by the LHC SUSY cross section working group at NLO including next-to-leading logarithms [63,64]. For a given model, new particles are excluded at 95% CL in mass regions where the theoretical prediction lies at or above the observed upper limit for the appropriate final state in Figs. 10-12. From Fig. 10, the 95% CL lower bound on the mass of a benchmark model of a diquark, for α true values between 0.11 and 0.42, varies between 7.9 and 6.7 TeV (observed), and between 7.9 and 6.4 TeV (expected), respectively. From Fig. 11, for α true = 0.25 we exclude at 95% CL diquark masses less than 7.6 TeV (observed) and 7.8 TeV (expected). These limits are for a benchmark choice of coupling values of the diquark to pairs of up quarks, y uu = 0.4, and to pairs of vector-like quarks, y χ = 0.6. Other choices of coupling values can yield different values for the diquark cross sections and the mass limits [22]. From Fig. 12, we exclude at 95% CL top squarks with masses between 0.50 and 0.52 TeV, and between 0.58 and 0.77 TeV (observed), and less than 0.67 TeV (expected).
The observed limits in the resonant and the nonresonant search are highly correlated, because they use the same events to search for localized excesses (bumps), though they search for different signals. Figure 12 shows that in the nonresonant search, at a dijet resonance mass of roughly 0.95 TeV, the observed limit is more than 2 standard deviations (s.d.) above the expected limit. This is the direct result of the pulls near an average dijet mass of 1 TeV for 0.25 < α < 0.35 in Fig. 9 which resulted in a low p-value for the background fit, as previously discussed. There is a correlated effect in the resonant search. In Fig. 10 we see that the observed limit is also more than 2 s.d. above the expected limit for a four-jet resonance mass between 3 and 4 TeV for α true = 0.27, 0.29, and 0.31.   For the nonresonant search, Fig. 13 shows the local p-value as a function of dijet resonance mass. The most significant signal hypothesis occurs at a dijet resonance mass of 0.95 TeV, for which the local significance is 3.6 s.d. and the global significance is 2.5 s.d.. Here the global significance takes into account the look-elsewhere effect (LEE) [54], using pseudo-experiments to calculate the probability of the background hypothesis to produce a signal-like effect with at least the observed local significance, anywhere within the sensitive range of dijet resonance mass. For the resonant search, Fig. 14 shows the local p-value as a function of four-jet resonance mass with α true = 0.25. The most significant signal hypothesis occurs at a four-jet resonance mass of 8.6 TeV and a dijet resonance mass of 2.1 TeV, for which the local and global significance are 3.9 s.d. and 1.6 s.d., respectively. Here again the global significance takes into account the LEE, using pseudo-experiments as described for the nonresonant search, but the LEE for the resonant search considers a significantly larger number of possible signals because of the 2D sensitive range of four-jet resonance mass and α true . At a four-jet resonance mass of 8.6 TeV the best fit value of σ B A is (1.56 +1. 40 −0.87 ) 10 −5 pb, which corresponds to 2.1 +2.0 −1.2 events. The second largest local significance is 2.9 s.d., for a four-jet resonance mass of ∼3 TeV and α true = 0.17. We note that the LEE assumes equal significance to an effect for any value of resonance mass, and the resonance width being smaller in absolute units at low masses yields higher probability for an effect within the low mass region.

Candidate Events
The high local significance, and corresponding global significance, of the effect observed at 8.6 TeV is predominantly caused by the two events which we previously identified, and we discuss them in more detail here and in the next two paragraphs. They are in Fig. 2, with a four-jet mass of roughly 8 TeV and an average dijet mass of roughly 2 TeV, and in Fig. 3 with the same four-jet mass. To better visualize the comparison of this small number of events to the distributions of the rest of the events, the background, and the signal hypothesis, we show in Fig. 15 the inclusive four-jet mass distribution corresponding to the sum of all thirteen α bins. This sums the two events which appear in Fig. 7 in two separate plots corresponding to α bins with 0.22 < α < 0.24 and 0.26 < α < 0.28. In this case, the backgroundonly fits are performed using the same functional forms, but with five parameters due to the larger event samples. These are the modified dijet function (ModDijet-5p) given by dσ/dm = p 0 (1 − x 1/3 ) p 1 /(x p 2 +p 3 log x+p 4 log 2 x ), the dijet function (Dijet-5p) given by dσ/dm = p 0 (1 − x) p 1 /(x p 2 +p 3 log x+p 4 log 2 x ), and the power-law times exponential function (PowExp-5p) given by dσ/dm = p 0 e −p 1 x−p 2 x 2 −p 3 x 3 /x p 4 , where x = m/ √ s, with m being the four-jet mass.
The first event, shown in Fig. 16, was collected during 2017, and has a four-jet mass of 8.0 TeV, while the invariant mass of both dijet pairs is 1.9 TeV. This is the event found in the CMS search in the dijet final state [1] which provided some of the motivation for this search. The p T of the first dijet pair, composed of jets labeled 1 and 3, is 3.5 TeV, and the second pair, consisting of jets 2 and 4, also has a p T of 3.5 TeV. The two dijet pairs are back-to-back in azimuthal angle (|∆φ| = 3.1), nearby in pseudorapidity (|∆η| = 0.4) and have a small mass asymmetry of 0.005. The two jets in each pair have a small distance with each other in the η-φ plane (∆R 1 = 1.0, ∆R 2 = 1.0), and the p T , η, and φ of each jet are shown in Fig. 16. Since the previous search [1], a recalibration of the jet energy scale caused the p T and mass values presented for this event to increase by roughly 1%.
The second event, shown in Fig. 17, was collected during 2018, and has a four-jet mass of 7.9 TeV, while the invariant mass of the first dijet pair that includes jets 1 and 4 is 2.0 TeV, and the mass of the second pair from jets 2 and 3 is 2.1 TeV. The two dijet pairs are back-to-back in azimuthal angle (|∆φ| = 3.1), nearby in pseudorapidity (|∆η| = 1.06), and have a mass asymmetry of 0.02. The ∆R 1 and ∆R 2 values between the two jets in the two dijet pairs are 1.5 and 1.3. Again, the p T , η, and φ of each jet can be seen in Fig. 17.

Summary
A search for paired dijet resonances has been presented in final states with at least four jets, probing dijet masses above 0.35 TeV and four-jet masses above 1.6 TeV. Data from protonproton collisions at √ s = 13 TeV were used in this search, collected by the CMS experiment at the LHC, corresponding to an integrated luminosity of 138 fb −1 .
For the first time, a search for resonant production of pairs of dijet resonances, where a massive intermediate state leads to a four-jet resonance in the final state, has been performed. A search for pairs of equal mass dijet resonances produced by a nonresonant process has also been conducted. Empirical functions that model the background, and simulated shapes of resonance signals, are fit to the observed four-jet and dijet mass distributions. There are two events in the tails of the distributions, with a four-jet mass of 8 TeV and an average dijet mass of 2 TeV. They result in a local significance of 3.9 standard deviations and a global significance of 1.6 standard deviations, when interpreted as a signal of a four-jet resonance with mass 8.6 TeV and a dijet resonance with mass 2.15 TeV. More data are needed to establish if these events are the first hints of a signal, and LHC data with the highest possible collision energy would be very effective.
Model-independent upper limits at 95% confidence level (CL) are presented on the production cross section times branching fraction and acceptance. For resonant production, the limits are presented as a function of the four-jet resonance mass between 2 and 9 TeV, for all accessible values of the ratio of the dijet to four-jet resonance masses. For nonresonant production, the limits are presented as functions of dijet resonance mass between 0.5 and 3.0 TeV.
Limits from the resonant search are compared to a model [22] of diquarks, which decay to pairs of vector-like quarks, which in-turn decay to a quark and a gluon. Mass limits for all accessible values of the ratio of the vector-like quark to diquark masses, for a benchmark scenario where the diquark couplings to pairs of up quarks is 0.4, and the diquark couplings to pairs of vectorlike quarks is 0.6 are presented. In particular, when the ratio of the vector-like quark to diquark mass is chosen to be 0.25, motivated by the two observed events discussed above, diquarks with a mass less than 7.6 TeV (observed) and 7.8 TeV (expected) are excluded at 95% CL.