1 Introduction

Particle jets with large transverse momenta \(p_{\mathrm {T}}\) are abundantly produced in proton–proton collisions at the CERN LHC through the strong interactions of quantum chromodynamics (QCD) between the incoming partons. When the momentum transfer is large, the dynamics can be predicted using perturbative techniques (pQCD). The two final-state partons at leading order (LO) in pQCD are produced back-to-back in the transverse plane, and thus the azimuthal angular separation between the two highest-\(p_{\mathrm {T}}\) jets, \(\varDelta \phi _\text {1,2} =|\phi _\text {jet1}-\phi _\text {jet2} |\), equals \(\pi \). The production of additional high-\(p_{\mathrm {T}}\) jets leads to a deviation of the azimuthal angle from \(\pi \). The measurement of azimuthal angular correlations (or decorrelation from \(\pi \)) in inclusive 2-jet topologies is a useful tool to test theoretical predictions of multijet production processes. Previous measurements of azimuthal correlation in inclusive 2-jet events were reported by the D0 Collaboration in \(\mathrm {p}\overline{\mathrm {p}} \) collisions at \(\sqrt{s}=1.96\,\text {Te}\text {V} \) at the Fermilab Tevatron [1, 2], and by the ATLAS Collaboration in \(\mathrm {p}\mathrm {p}\) collisions at \(\sqrt{s}=7\,\text {Te}\text {V} \) [3] and the CMS Collaboration in \(\mathrm {p}\mathrm {p}\) collisions at \(\sqrt{s}=7\) and \(8\,\text {Te}\text {V} \) [4, 5] at the LHC. Multijet correlations have been measured by the ATLAS Collaboration at \(\sqrt{s}=8\,\text {Te}\text {V} \) [6, 7].

This paper reports measurements of the normalized inclusive 2-, 3-, and 4-jet cross sections as a function of the azimuthal angular separation between the two highest \(p_{\mathrm {T}}\) (leading) jets, \(\varDelta \phi _\text {1,2}\),

$$\begin{aligned} \frac{1}{\sigma } \frac{\mathrm {d}\sigma }{\mathrm {d}\varDelta \phi _\text {1,2}}, \end{aligned}$$

for several regions of the leading jet \(p_{\mathrm {T}}\), \(p_{\mathrm {T}} ^{\text {max}}\), for the rapidity region \(|y |<2.5\). The measurements cover the region \(\pi /2 < \varDelta \phi _\text {1,2} \le \pi \); the region \(\varDelta \phi _\text {1,2} \le \pi /2\) includes large backgrounds due to \({\mathrm {t}\overline{\mathrm {t}}}\) and Z/W+jet(s) events. Experimental and theoretical uncertainties are reduced by normalizing the \(\varDelta \phi _\text {1,2}\) distribution to the total dijet cross section within each region of \(p_{\mathrm {T}} ^{\text {max}}\).

For 3- and 4-jet topologies, measurements of the normalized inclusive 3- and 4-jet cross sections are also presented as a function of the minimum azimuthal angular separation between any two of the three or four highest \(p_{\mathrm {T}}\) jets, \(\varDelta \phi _\text {2j}^\text {min}\),

$$\begin{aligned} \frac{1}{\sigma } \frac{\mathrm {d}\sigma }{\mathrm {d}\varDelta \phi _\text {2j}^\text {min}}, \end{aligned}$$

for several regions of \(p_{\mathrm {T}} ^{\text {max}}\), for \(|y |<2.5\). This observable, which is infrared safe (independent of additional soft radiation), is especially suited for studying correlations amongst the jets in multijet events: the maximum value of \(\varDelta \phi _\text {2j}^\text {min}\) is \(2 \pi /3\) for 3-jet events (the “Mercedes star” configuration), while it is \(\pi /2\) in the 4-jet case (corresponding to the “cross” configuration). The cross section for small angular separations is suppressed because of the finite jet sizes for a particular jet algorithm. The observable \(\varDelta \phi _\text {2j}^\text {min}\) is sensitive to the contributions of jets with lower \(p_{\mathrm {T}}\) than the leading jet, i.e. the subleading jets, and one can distinguish nearby (nearly collinear) jets (at large \(\varDelta \phi _\text {2j}^\text {min}\)) from other additional high \(p_{\mathrm {T}}\) jets (small \(\varDelta \phi _\text {2j}^\text {min}\)), yielding information additional to that of the \(\varDelta \phi _\text {1,2}\) observable. The 4-jet cross section differential in \(\varDelta \phi _\text {2j}^\text {min}\) has also been measured by the ATLAS Collaboration [7].

The measurements are performed using data collected during 2016 with the CMS experiment at the LHC, and the event sample corresponds to an integrated luminosity of 35.9\(\,\text {fb}^{-1}\) of proton–proton collisions at \(\sqrt{s}=13\,\text {Te}\text {V} \).

2 The CMS detector

The central feature of the CMS detector is a superconducting solenoid, \(13~\hbox {m}\) in length and \(6~\hbox {m}\) in inner diameter, providing an axial magnetic field of \(3.8~\hbox {T}\). Within the solenoid volume are a silicon pixel and strip tracker, a lead tungstate crystal electromagnetic calorimeter (ECAL) and a brass and scintillator hadron calorimeter (HCAL), each composed of a barrel and two endcap sections. Charged-particle trajectories are measured by the tracker with full azimuthal coverage within pseudorapidities \(|\eta |< 2.5\). The ECAL, which is equipped with a preshower detector in the endcaps, and the HCAL cover the region \(|\eta |< 3.0\). Forward calorimeters extend the pseudorapidity coverage provided by the barrel and endcap detectors to the region \(3.0< |\eta | < 5.2\). Finally, muons are measured up to \(|\eta |< 2.4\) by gas-ionization detectors embedded in the steel flux-return yoke outside the solenoid. A 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. [8].

3 Theoretical predictions

Predictions from five different Monte Carlo (MC) event generators are compared with data. The pythia  8  [9] and herwig++  [10] event generators are used, both based on LO \(2\rightarrow 2\) matrix element calculations. The pythia  8 event generator simulates parton showers ordered in \(p_{\mathrm {T}}\) and uses the Lund string model [11] for hadronization, while herwig++ generates parton showers through angular-ordered emissions and uses a cluster fragmentation model [12] for hadronization. The contribution of multiparton interactions (MPI) is simulated in both pythia  8 and herwig++, but the number of generated MPI varies between pythia  8 and herwig++ MPI simulations. The MPI parameters of both generators are tuned to measurements in proton–proton collisions at the LHC and proton–antiproton collisions at the Tevatron [13], while the hadronization parameters are determined from fits to LEP data. For pythia  8 the CUETP8M1 [13] tune, which is based on the NNPDF2.3LO PDF set [14, 15], is employed, while for herwig++ the CUETHppS1 tune [13], based on the CTEQ6L1 PDF set [16], is used.

The MadGraph  [17, 18] event generator provides LO matrix element calculations with up to four outgoing partons, i.e. \(2\rightarrow 2\), \(2\rightarrow 3\), and \(2\rightarrow 4\) diagrams. It is interfaced to pythia  8 with tune CUETP8M1 for the implementation of parton showers, hadronization, and MPI. In order to match with pythia  8 the \(k_{\mathrm {T}}\)-MLM matching procedure [19] with a matching scale of \(14 \,\text {Ge}\text {V} \) is used to avoid any double counting of the parton configurations generated within the matrix element calculation and the ones simulated by the parton shower. The NNPDF2.3LO PDF set is used for the hard-process calculation.

Predictions based on next-to-leading-order (NLO) pQCD are obtained with the powheg box library [20,21,22] and the herwig  7  [23] event generator. The events simulated with powheg are matched to pythia  8 or to herwig++ parton showers and MPI, while herwig  7 uses similar parton shower and MPI models as herwig++, and the MC@NLO [24, 25] method is applied to combine the parton shower with the NLO calculation. The powheg generator is used in the NLO dijet mode [26], referred to as ph-2j, as well as in the NLO three-jet mode [27], referred to as ph-3j, both using the NNPDF3.0NLO PDF set [28]. The powheg generator, referred to as ph-2j-lhe, is also used in the NLO dijet mode without parton showers and MPI. A minimum \(p_{\mathrm {T}}\) for real parton emission of 10 GeV is required for the ph-2j predictions, and similarly for the ph-3j predictions a minimum \(p_{\mathrm {T}}\) for the three final-state partons of 10 \(\,\text {Ge}\text {V}\) is imposed. To simulate the contributions due to parton showers, hadronization, and MPIs, the ph-2j is matched to pythia  8 with tune CUETP8M1 and herwig++ with tune CUETHppS1, while the ph-3j is matched only to pythia  8 with tune CUETP8M1. The matching between the powheg matrix element calculations and the pythia  8 underlying event (UE) simulation is performed using the shower-veto procedure, which rejects showers if their transverse momentum is greater than the minimal \(p_{\mathrm {T}}\) of all final-state partons simulated in the matrix element (parameter pthard = 2 [26]). Predictions from the herwig  7 event generator are based on the MMHT2014 PDF set [29] and the default tune H7-UE-MMHT [23] for the UE simulation. A summary of the details of the MC event generators used for comparisons with the experimental data is shown in Table 1.

Table 1 Monte Carlo event generators used for comparison in this analysis. Version of the generators, PDF set, underlying event tune, and corresponding references are listed

Uncertainties in the theoretical predictions of the parton shower simulation are illustrated using the pythia  8 event generator. Choices of scale for the parton shower are expected to have the largest impact on the azimuthal distributions. The parton shower uncertainty is calculated by independently varying the renormalization scales (\(\mu _r\)) for initial- and final-state radiation by a factor 2 in units of the \(p_{\mathrm {T}}\) of the emitted partons of the hard scattering. The maximum deviation found is considered a theoretical uncertainty in the event generator predictions.

4 Jet reconstruction and event selection

The measurements are based on data samples collected with single-jet high-level triggers (HLT) [30, 31]. Five such triggers are considered that require at least one jet in an event with \(p_{\mathrm {T}} > 140\), 200, 320, 400, or \(450\,\text {Ge}\text {V} \) in the full rapidity coverage of the CMS detector. All triggers are prescaled except the one with the highest threshold. Table 2 shows the integrated luminosity \(\mathcal {L}\) for the five trigger samples. The relative efficiency of each trigger is estimated using triggers with lower \(p_{\mathrm {T}}\) thresholds. Using these five jet energy thresholds, a 100% trigger efficiency is achieved in the region of \(p_{\mathrm {T}} ^{\text {max}} > 200\,\text {Ge}\text {V} \).

Table 2 The integrated luminosity for each trigger sample considered in this analysis

Particles are reconstructed and identified using a particle-flow (PF) algorithm [32], which uses an optimized combination of information from the various elements of the CMS detector. Jets are reconstructed by clustering the Lorentz vectors of the PF candidates with the infrared- and collinear-safe anti-\(k_{\mathrm {T}}\) clustering algorithm [33] with a distance parameter \(R=0.4\). The clustering is performed with the FastJet package [34]. The technique of charged-hadron subtraction [35] is used to remove tracks identified as originating from additional pp interactions within the same or neighbouring bunch crossings (pileup). The average number of pileup interactions observed in the data is about 27.

The reconstructed jets require energy corrections to account for residual nonuniformities and nonlinearities in the detector response. These jet energy scale (JES) corrections [35] are derived using simulated events that are generated with pythia  8.219 [9] using tune CUETP8M1 [13] and processed through the CMS detector simulation based on Geant4 [36]; they are confirmed with in situ measurements with dijet, multijet, photon+jet, and leptonic Z+jet events. An offset correction is required to account for the extra energy clustered into jets due to pileup. The JES corrections, which depend on the \(\eta \) and \(p_{\mathrm {T}}\) of the jet, are applied as multiplicative factors to the jet four-momentum vectors. The typical overall correction is about 10% for central jets having \(p_{\mathrm {T}} = 100\,\text {Ge}\text {V} \) and decreases with increasing \(p_{\mathrm {T}}\).

Resolution studies on the measurements of \(\varDelta \phi _\text {1,2}\) and \(\varDelta \phi _\text {2j}^\text {min}\) are performed using pythia  8.219 with tune CUETP8M1 processed through the CMS detector simulation. The azimuthal angular separation is determined with an accuracy from \(1^\circ \) to \(0.5^\circ \) (0.017 to 0.0087 in radians) for \(p_{\mathrm {T}} ^{\text {max}} = 200\,\text {Ge}\text {V} \) to \(1\,\text {Te}\text {V} \), respectively.

Events are required to have at least one primary vertex candidate [37] reconstructed offline from at least five charged-particle tracks and lies along the beam line within 24\(\text {\,cm}\) of the nominal interaction point. The reconstructed vertex with the largest value of summed physics-object \(p_{\mathrm {T}} ^2\) is taken to be the primary \(\mathrm {p}\mathrm {p}\) interaction vertex. The physics objects are the objects determined by a jet finding algorithm [33, 34] applied to all charged tracks associated with the vertex plus the corresponding associated missing transverse momentum. Additional selection criteria are applied to each event to remove spurious jet-like signatures originating from isolated noise patterns in certain HCAL regions. Stringent criteria [38] are applied to suppress these nonphysical signatures; each jet should contain at least two particles, one of which is a charged hadron, and the jet energy fraction carried by neutral hadrons and photons should be less than 90%. These criteria have a jet selection efficiency greater than 99% for genuine jets.

For the measurements of the normalized inclusive 2-, 3-, and 4-jet cross sections as a function of \(\varDelta \phi _\text {1,2}\) or \(\varDelta \phi _\text {2j}^\text {min}\) all jets in the event with \(p_{\mathrm {T}} > 100 \,\text {Ge}\text {V} \) and a rapidity \(|y |<5\) are considered and ordered in \(p_{\mathrm {T}}\). Events are selected where the two highest-\(p_{\mathrm {T}}\) jets have \(|y |<2.5\), (i.e. events are not counted where one of the leading jets has \(|y |>2.5\)). Also, events are only selected in which the highest-\(p_{\mathrm {T}}\) jet has \(|y |<2.5\) and exceeds 200\(\,\text {Ge}\text {V}\). The inclusive 2-jet event sample includes events where the two leading jets lie within the tracker coverage of \(|y |<2.5\). Similarly the 3-jet (4-jet) event sample includes those events where the three (four) leading jets lie within \(|y |<2.5\), respectively. In this paper results are presented in bins of \(p_{\mathrm {T}} ^{\text {max}}\), corresponding to the \(p_{\mathrm {T}}\) of the leading jet, which is always within \(|y |<2.5\).

5 Measurements of the normalized inclusive 2-, 3-, and 4-jet cross sections in \(\varDelta \phi _\text {1,2}\) and \(\varDelta \phi _\text {2j}^\text {min}\)

The normalized inclusive 2-, 3-, and 4-jet cross sections differential in \(\varDelta \phi _\text {1,2}\) and \(\varDelta \phi _\text {2j}^\text {min}\) are corrected for the finite detector resolution to better approximate the final-state particles, a procedure called “unfolding”. In this way, a direct comparison of this measurement to results from other experiments and to QCD predictions is possible. Particles are considered stable if their mean decay length is \(c\tau > 1\text {\,cm} \).

The bin width used in the measurements of \(\varDelta \phi _\text {1,2}\) and \(\varDelta \phi _\text {2j}^\text {min}\) is set to \(\pi / 36 = 0.087\) rads (\(5^\circ \)), which is five to ten times larger than the azimuthal angular separation resolution. The corrections due to the unfolding are approximately a few per cent.

The unfolding procedure is based on the matrix inversion algorithm implemented in the software package RooUnfold [39] using a 2-dimensional response matrix that correlates the modeled distribution with the reconstructed one. The response matrix is created by the convolution of the \(\varDelta \phi \) resolution with the generator-level inclusive 2-, 3-, and 4- cross section distributions from pythia  8 with tune CUETP8M1. The unfolded distributions differ from the distributions at detector level by 1–4%. As a cross-check, the above procedure was repeated by creating the response matrix with event samples obtained with the full Geant4 detector simulation, and no significant difference was observed.

We consider three main sources of systematic uncertainties that arise from the estimation of the JES calibration, the jet energy resolution (JER), and the unfolding correction. The relative JES uncertainty is estimated to be 1–2% for PF jets using charged-hadron subtraction [35]. The resulting uncertainties in the normalized 2-, 3-, and 4-jet cross sections differential in \(\varDelta \phi _\text {1,2}\) range from 3% at \(\pi /2\) to 0.1% at \(\pi \). For the normalized 3- and 4-jet cross sections differential in \(\varDelta \phi _\text {2j}^\text {min}\) the resulting uncertainties range from 0.1 to 1%, and 0.1–2%, respectively.

The JER [35] is responsible for migration of events among the \(p_{\mathrm {T}} ^{\text {max}}\) regions, and its parametrization is determined from a full detector simulation using events generated by pythia  8 with tune CUETP8M1. The effect of the JER uncertainty is estimated by varying its parameters within their uncertainties [35] and comparing the normalized inclusive 2-, 3-, and 4-jet cross sections before and after the changes. The JER-induced uncertainty ranges from 1% at \(\pi /2\) to 0.1% at \(\pi \) for the normalized 2-, 3-, and 4-jet cross sections differential in \(\varDelta \phi _\text {1,2}\) and is less than 0.5% for the normalized 3- and 4-jet cross sections differential in \(\varDelta \phi _\text {2j}^\text {min}\).

The above systematic uncertainties in the JES calibration and the JER cover the effects from migrations due to the \(p_{\mathrm {T}}\) thresholds, i.e. migrations between the 2-, 3-, and 4-jet samples and migrations between the various \(p_{\mathrm {T}} ^{\text {max}}\) regions of the measurements.

The unfolding procedure is affected by uncertainties in the parametrization of the \(\varDelta \phi \) resolution. Alternative response matrices, generated by varying the \(\varDelta \phi \) resolution by ±10%, are used to unfold the measured spectra. This variation is motivated by studies on the \(\varDelta \phi \) resolution for simulated di-jet events [32]. The uncertainty in the unfolding correction factors is estimated to be about 0.2%. An additional systematic uncertainty is obtained by examining the dependence of the response matrix on the choice of the MC generator. Alternative response matrices are constructed using the herwig++ event generator [10] with tune EE5C [40]; the effect is <0.1%. A total systematic unfolding uncertainty of 0.2% is considered, which accounts for all these various uncertainty sources.

6 Comparison with theoretical predictions

6.1 The \(\varDelta \phi _\text {1,2}\) measurements

The unfolded, normalized, inclusive 2-, 3-, and 4-jet cross sections differential in \(\varDelta \phi _\text {1,2}\) are shown in Figs. 1, 2, 3 for the various \(p_{\mathrm {T}} ^{\text {max}}\) regions considered in this analysis. In the 2-jet case the \(\varDelta \phi _\text {1,2}\) distributions are strongly peaked at \(\pi \) and become steeper with increasing \(p_{\mathrm {T}} ^{\text {max}}\). In the 3-jet case, the \(\varDelta \phi _\text {1,2}\) distributions become flatter at \(\pi \), since by definition dijet events do not contribute, and in the 4-jet case they become even flatter. The data points are overlaid with the predictions from the ph-2j + pythia  8 event generator.

Fig. 1
figure 1

Normalized inclusive 2-jet cross section differential in \(\varDelta \phi _\text {1,2}\) for nine \(p_{\mathrm {T}} ^{\text {max}}\) regions, scaled by multiplicative factors for presentation purposes. The size of the data symbol includes both statistical and systematic uncertainties. The data points are overlaid with the predictions from the ph-2j + pythia  8 event generator

Fig. 2
figure 2

Normalized inclusive 3-jet cross section differential in \(\varDelta \phi _\text {1,2}\) for eight \(p_{\mathrm {T}} ^{\text {max}}\) regions, scaled by multiplicative factors for presentation purposes. The size of the data symbol includes both statistical and systematic uncertainties. The data points are overlaid with the predictions from the ph-2j + pythia  8 event generator

Fig. 3
figure 3

Normalized inclusive 4-jet cross section differential in \(\varDelta \phi _\text {1,2}\) for eight \(p_{\mathrm {T}} ^{\text {max}}\) regions, scaled by multiplicative factors for presentation purposes. The size of the data symbol includes both statistical and systematic uncertainties. The data points are overlaid with the predictions from the ph-2j + pythia  8 event generator

The ratios of the pythia  8, herwig++, and MadGraph + pythia  8 event generator predictions to the normalized inclusive 2-, 3-, and 4-jet cross section differential in \(\varDelta \phi _\text {1,2}\) are shown in Figs. 4, 5, and 6, respectively, for all \(p_{\mathrm {T}} ^{\text {max}}\) regions. The solid band around unity represents the total experimental uncertainty and the error bars on the points represent the statistical uncertainties in the simulated data. Among the LO dijet event generators, herwig++ exhibits the largest deviations from the experimental measurements, whereas pythia  8 behaves much better than herwig++, although with deviations of up to 30–40%, in particular around \( \varDelta \phi _\text {1,2} = 5\pi /6\) in the 2-jet case and around \(\varDelta \phi _\text {1,2} < 2\pi /3\) in the 3- and 4-jet case. Predictions from herwig++ tend to overestimate the measurements as a function of \( \varDelta \phi _\text {1,2} \) in the 2-, 3-, and 4-jet cases, especially at \( \varDelta \phi _\text {1,2} < 5\pi /6\) for \(p_{\mathrm {T}} ^{\text {max}} > 400 \,\text {Ge}\text {V} \). However, it is remarkable that predictions based on the \(2\rightarrow 2\) matrix element calculations supplemented with parton showers, MPI, and hadronization describe the \(\varDelta \phi _\text {1,2} \) distributions rather well, even in regions that are sensitive to hard jets not included in the matrix element calculations. The MadGraph + pythia  8 calculation using up to 4 partons in the matrix element calculations provides the best description of the measurements.

Fig. 4
figure 4

Ratios of pythia  8, herwig++, and MadGraph + pythia  8 predictions to the normalized inclusive 2-jet cross section differential in \(\varDelta \phi _\text {1,2}\), for all \(p_{\mathrm {T}} ^{\text {max}}\) regions. The solid band indicates the total experimental uncertainty and the vertical bars on the points represent the statistical uncertainties in the simulated data

Fig. 5
figure 5

Ratios of pythia  8, herwig++, and MadGraph + pythia  8 predictions to the normalized inclusive 3-jet cross section differential in \(\varDelta \phi _\text {1,2}\), for all \(p_{\mathrm {T}} ^{\text {max}}\) regions. The solid band indicates the total experimental uncertainty and the vertical bars on the points represent the statistical uncertainties in the simulated data

Fig. 6
figure 6

Ratios of pythia  8, herwig++, and MadGraph + pythia  8 predictions to the normalized inclusive 4-jet cross section differential in \(\varDelta \phi _\text {1,2}\), for all \(p_{\mathrm {T}} ^{\text {max}}\) regions. The solid band indicates the total experimental uncertainty and the vertical bars on the points represent the statistical uncertainties in the simulated data

Figures 7, 8 and 9 show the ratios of the ph-2j matched to pythia  8 and herwig++, ph-3j + pythia  8, and herwig  7 event generators predictions to the normalized inclusive 2-, 3-, and 4-jet cross section differential in \(\varDelta \phi _\text {1,2}\), for all \(p_{\mathrm {T}} ^{\text {max}}\) regions. The solid band around unity represents the total experimental uncertainty and the vertical bars on the points represent the statistical uncertainties in the simulated data. The predictions of ph-2j and ph-3j exhibit deviations from the measurement, increasing towards small \(\varDelta \phi _\text {1,2}\). While ph-2j is above the data, ph-3j predicts too few events at small \(\varDelta \phi _\text {1,2}\). These deviations were investigated in a dedicated study with parton showers and MPI switched off. Because of the kinematic restriction of a 3-parton state, ph-2j without parton showers cannot fill the region \(\varDelta \phi _\text {1,2} <2\pi /3\), shown as ph-2j-lhe with the dashed line in Fig. 7, whereas for ph-3j the parton showers have little impact. Thus, the events at low \(\varDelta \phi _\text {1,2}\) observed for ph-2j originate from leading-log parton showers, and there are too many of these. In contrast, the ph-3j prediction, which provides \(2\rightarrow 3\) jet calculations at NLO QCD, is below the measurement. The NLO ph-2j calculation and the LO powheg three-jet calculation are equivalent when initial- and final-state radiation are not allowed to occur.

Fig. 7
figure 7

Ratios of ph-2j + pythia  8, ph-2j-lhe, ph-2j + herwig++, ph-3j + pythia  8, and herwig  7 predictions to the normalized inclusive 2-jet cross section differential in \(\varDelta \phi _\text {1,2}\), for all \(p_{\mathrm {T}} ^{\text {max}}\) regions. The solid band indicates the total experimental uncertainty and the vertical bars on the points represent the statistical uncertainties in the simulated data

Fig. 8
figure 8

Ratios of ph-2j + pythia  8, ph-2j + herwig++, ph-3j + pythia  8, and herwig  7 predictions to the normalized inclusive 3-jet cross section differential in \(\varDelta \phi _\text {1,2}\), for all \(p_{\mathrm {T}} ^{\text {max}}\) regions. The solid band indicates the total experimental uncertainty and the vertical bars on the points represent the statistical uncertainties in the simulated data

Fig. 9
figure 9

Ratios of ph-2j + pythia  8, ph-2j + herwig++, ph-3j + pythia  8, and herwig  7 predictions to the normalized inclusive 4-jet cross section differential in \(\varDelta \phi _\text {1,2}\), for all \(p_{\mathrm {T}} ^{\text {max}}\) regions. The solid band indicates the total experimental uncertainty and the vertical bars on the points represent the statistical uncertainties in the simulated data

The predictions from ph-2j matched to pythia  8 describe the normalized cross sections better than those where ph-2j is matched to herwig++. Since the hard process calculation is the same, the difference between the two predictions might be due to the treatment of parton showers in pythia  8 and herwig++ and to the matching to the matrix element calculation. The pythia  8 and herwig++ parton shower calculations use different \(\alpha _S\) values for initial- and final-state emissions, in addition to a different upper scale for the parton shower simulation, which is higher in pythia  8 than in herwig++. The dijet NLO calculation of herwig  7 provides the best description of the measurements, indicating that the MC@NLO method of combining parton showers with the NLO parton level calculations has advantages compared to the POWHEG method in this context.

For \(\varDelta \phi _\text {1,2} \) generator-level predictions in the 2-jet case, parton shower uncertainties have a very small impact (<5%) at values close to \(\pi \) and go up to 40–60% for increasing \(p_{\mathrm {T}} ^{\text {max}} \) at \(\varDelta \phi _\text {1,2} \sim \pi /2\). For the 3- and 4-jet scenarios, parton shower uncertainties are less relevant, not exceeding \(\sim \)20% for \(\varDelta \phi _\text {1,2} \).

6.2 The \(\varDelta \phi _\text {2j}^\text {min}\) measurements

The unfolded, normalized, inclusive 3- and 4-jet cross sections differential in \(\varDelta \phi _\text {2j}^\text {min}\) are shown in Figs. 10 and 11, respectively, for eight \(p_{\mathrm {T}} ^{\text {max}}\) regions. The measured distributions decrease towards the kinematic limit of \(\varDelta \phi _\text {2j}^\text {min} \rightarrow 2 \pi /3( \pi /2)\) for the 3-jet and 4-jet case, respectively. The data points are overlaid with the predictions from the ph-2j + pythia  8 event generator. The size of the data symbol includes both statistical and systematic uncertainties.

Fig. 10
figure 10

Normalized inclusive 3-jet cross section differential in \(\varDelta \phi _\text {2j}^\text {min}\) for eight \(p_{\mathrm {T}} ^{\text {max}}\) regions, scaled by multiplicative factors for presentation purposes. The size of the data symbol includes both statistical and systematic uncertainties. The data points are overlaid with the predictions from the ph-2j + pythia  8 event generator

Fig. 11
figure 11

Normalized inclusive 4-jet cross section differential in \(\varDelta \phi _\text {2j}^\text {min}\) for eight \(p_{\mathrm {T}} ^{\text {max}}\) regions, scaled by multiplicative factors for presentation purposes. The size of the data symbol includes both statistical and systematic uncertainties. The data points are overlaid with the predictions from the ph-2j + pythia  8 event generator

Figures 12 and 13 show, respectively, the ratios of the pythia  8, herwig++, and MadGraph + pythia  8 event generators predictions to the normalized inclusive 3- and 4-jet cross sections differential in \(\varDelta \phi _\text {2j}^\text {min}\), for all \(p_{\mathrm {T}} ^{\text {max}}\) regions. The pythia  8 event generator shows larger deviations from the measured \(\varDelta \phi _\text {2j}^\text {min}\) distributions in comparison to herwig++, which provides a reasonable description of the measurement. The MadGraph generator matched to pythia  8 provides a reasonable description of the measurements in the 3-jet case, but shows deviations in the 4-jet case.

The predictions from MadGraph + pythia  8 and pythia  8 are very similar for the normalized cross sections as a function of \(\varDelta \phi _\text {2j}^\text {min}\) in the four-jet case. It has been checked that predictions obtained with the MadGraph matrix element with up to 4 partons included in the calculation without contribution of the parton shower are able to reproduce the data very well. Parton shower effects increase the number of events with low values of \(\varDelta \phi _\text {2j}^\text {min}\).

Fig. 12
figure 12

Ratios of pythia  8, herwig++, and MadGraph + pythia  8 predictions to the normalized inclusive 3-jet cross section differential in \(\varDelta \phi _\text {2j}^\text {min}\), for all \(p_{\mathrm {T}} ^{\text {max}}\) regions. The solid band indicates the total experimental uncertainty and the vertical bars on the points represent the statistical uncertainties in the simulated data

Fig. 13
figure 13

Ratios of pythia  8, herwig++, and MadGraph + pythia  8 predictions to the normalized inclusive 4-jet cross section differential in \(\varDelta \phi _\text {2j}^\text {min}\), for all \(p_{\mathrm {T}} ^{\text {max}}\) regions. The solid band indicates the total experimental uncertainty and the vertical bars on the points represent the statistical uncertainties in the simulated data

Figures 14 and 15 illustrate the ratios of predictions from ph-2j matched to pythia  8 and herwig++, ph-3j + pythia  8, and herwig  7 to the normalized inclusive 3- and 4-jet cross sections differential in \(\varDelta \phi _\text {2j}^\text {min}\), for all \(p_{\mathrm {T}} ^{\text {max}}\) regions. Due to an unphysical behavior of the herwig  7 prediction (which has been confirmed by the herwig  7 authors), the first \(\varDelta \phi _\text {2j}^\text {min}\) and last \(\varDelta \phi _\text {1,2}\) bins are not shown in Figs. 8, 9, 14 and 15. An additional uncertainty is introduced to the prediction of herwig  7, that is evaluated as the difference between this prediction and the prediction when the first bin is replaced with the result from herwig++. The additional uncertainty ranges from 2 to 10%. Among the three NLO dijet calculations ph-2j matched to pythia  8 or to herwig++ provides the best description of the measurements.

For the two lowest \(p_{\mathrm {T}} ^{\text {max}}\) regions in Figs. 13 and 15, which correspond to the 4-jet case, the measurements become statistically limited because the data used for these two regions were collected with highly prescaled triggers with \(p_{\mathrm {T}}\) thresholds of 140 and 200 \(\,\text {Ge}\text {V}\) (c.f. Table 2).

The ph-3j predictions suffer from low statistical accuracy, especially in the highest interval of \(p_{\mathrm {T}} ^{\text {max}}\), because the same \(p_{\mathrm {T}}\) threshold is applied to all 3 jets resulting in low efficiency at large \(p_{\mathrm {T}}\). Nevertheless, the performance of the ph-3j simulation on multijet observables can already be inferred by the presented predictions, especially in the low \(p_{\mathrm {T}}\) region.

The effect of parton shower uncertainties in the event generator predictions of \(\varDelta \phi _\text {2j}^\text {min} \) is estimated to be less than 10% over the entire phase space.

Fig. 14
figure 14

Ratios of ph-2j + pythia  8, ph-2j + herwig++, ph-3j + pythia  8, and herwig  7 predictions to the normalized inclusive 3-jet cross section differential in \(\varDelta \phi _\text {2j}^\text {min}\), for all \(p_{\mathrm {T}} ^{\text {max}}\) regions. The solid band indicates the total experimental uncertainty and the vertical bars on the points represent the statistical uncertainties of the simulated data

Fig. 15
figure 15

Ratios of ph-2j + pythia  8, ph-2j + herwig++, ph-3j + pythia  8, and herwig  7 predictions to the normalized inclusive 4-jet cross section differential in \(\varDelta \phi _\text {2j}^\text {min}\), for all \(p_{\mathrm {T}} ^{\text {max}}\) regions. The solid band indicates the total experimental uncertainty and the vertical bars on the points represent the statistical uncertainties of the simulated data

7 Summary

Measurements of the normalized inclusive 2-, 3-, and 4-jet cross sections differential in the azimuthal angular separation \(\varDelta \phi _\text {1,2}\) and of the normalized inclusive 3- and 4-jet cross sections differential in the minimum azimuthal angular separation between any two jets \(\varDelta \phi _\text {2j}^\text {min}\) are presented for several regions of the leading-jet transverse momentum \(p_{\mathrm {T}} ^{\text {max}}\). The measurements are performed using data collected during 2016 with the CMS detector at the CERN LHC corresponding to an integrated luminosity of 35.9\(\,\text {fb}^{-1}\) of proton–proton collisions at \(\sqrt{s}=13\,\text {Te}\text {V} \).

The measured distributions in \(\varDelta \phi _\text {1,2}\) and \(\varDelta \phi _\text {2j}^\text {min}\) are compared with predictions from pythia  8, herwig++, MadGraph + pythia  8, ph-2j matched to pythia  8 and herwig++, ph-3j + pythia  8, and herwig  7 event generators.

The leading order (LO) pythia  8 dijet event generator exhibits small deviations from the \(\varDelta \phi _\text {1,2}\) measurements but shows significant deviations at low-\(p_{\mathrm {T}}\) in the \(\varDelta \phi _\text {2j}^\text {min}\) distributions. The herwig++ event generator exhibits the largest deviations of any of the generators for the \(\varDelta \phi _\text {1,2}\) measurements, but provides a reasonable description of the \(\varDelta \phi _\text {2j}^\text {min}\) distributions. The tree-level multijet event generator MadGraph in combination with pythia  8 for showering, hadronization, and multiparton interactions provides a good overall description of the measurements, except for the \(\varDelta \phi _\text {2j}^\text {min}\) distributions in the 4-jet case, where the generator deviates from the measurement mainly at high \(p_{\mathrm {T}} ^{\text {max}}\).

The dijet next-to-leading order (NLO) ph-2j event generator deviates from the \(\varDelta \phi _\text {1,2}\) measurements, but provides a good description of the \(\varDelta \phi _\text {2j}^\text {min}\) observable. The predictions from the three-jet NLO ph-3j event generator exhibit large deviations from the measurements and describe the considered multijet observables in a less accurate way than the predictions from ph-2j. Parton shower contributions are responsible for the different behaviour of the ph-2j and ph-3j predictions. Finally, predictions from the dijet NLO herwig  7 event generator matched to parton shower contributions with the MC@NLO method provide a very good description of the \(\varDelta \phi _\text {1,2}\) measurements, showing improvement in comparison to herwig++.

All these observations emphasize the need to improve predictions for multijet production. Similar observations, for the inclusive 2-jet cross sections differential in \(\varDelta \phi _\text {1,2}\), were reported previously by CMS [5] at a different centre-of-mass energy of 8\(\,\text {Te}\text {V}\). The extension of \(\varDelta \phi _\text {1,2}\) correlations, and the measurement of the \(\varDelta \phi _\text {2j}^\text {min}\) distributions in inclusive 3- and 4-jet topologies are novel measurements of the present analysis.