Determination of the strong coupling constant from transverse energy$-$energy correlations in multijet events at $\sqrt{s} = 13$ TeV with the ATLAS detector

Measurements of transverse energy$-$energy correlations and their associated azimuthal asymmetries in multijet events are presented. The analysis is performed using a data sample corresponding to 139 $\mbox{fb\(^{-1}\)}$ of proton$-$proton collisions at a centre-of-mass energy of $\sqrt{s} = 13$ TeV, collected with the ATLAS detector at the Large Hadron Collider. The measurements are presented in bins of the scalar sum of the transverse momenta of the two leading jets and unfolded to particle level. They are then compared to next-to-next-to-leading-order perturbative QCD calculations for the first time, which feature a significant reduction in the theoretical uncertainties estimated using variations of the renormalisation and factorisation scales. The agreement between data and theory is good, thus providing a precision test of QCD at large momentum transfers $Q$. The strong coupling constant $\alpha_s$ is extracted differentially as a function of $Q$, showing a good agreement with the renormalisation group equation and with previous analyses. A simultaneous fit to all transverse energy$-$energy correlation distributions across different kinematic regions yields a value of $\alpha_\mathrm{s}(m_Z) = 0.1175 \pm 0.0006 \mbox{ (exp.)} ^{+0.0034}_{-0.0017} \mbox{ (theo.)}$, while the global fit to the asymmetry distributions yields $\alpha_{\mathrm{s}}(m_Z) = 0.1185 \pm 0.0009 \mbox{ (exp.)} ^{+0.0025}_{-0.0012} \mbox{ (theo.)}$.


Introduction
Multĳet final states, produced in proton-proton ( ) collisions with large momentum transfer at the LHC provide an ideal testing ground for perturbative Quantum Chromodynamics (pQCD). Event shapes [1,2] are a class of observables defined as functions of the final-state particles four-momenta, which characterise the hadronic energy flow in a collision. They can be used to precisely test pQCD calculations and additionally to extract the value of the strong coupling constant, s . Event shape variables were measured in + − collision experiments from PETRA and PEP [3][4][5] to LEP and SLC [6][7][8][9][10], at the collider HERA [11][12][13][14][15] and in hadron-hadron collisions at the Tevatron [16] and the LHC [17][18][19].
Transverse energy-energy correlations (TEEC) and their associated azimuthal asymmetries (ATEEC) were proposed as the appropriate generalisation for hadron collider experiments in Ref. [41], where leading-order (LO) predictions were also presented. In experiments with incoming hadrons, observables that are invariant with respect to boosts along the direction of the beam axis are beneficial. As jet-based observables, the TEEC and ATEEC make use of the jet transverse energy T = /cosh , where is the jet energy and is the jet rapidity. The TEEC function is defined as the transverse-energy-weighted distribution of the azimuthal differences between jet pairs in the final state [42], i.e. 1 dΣ d cos where the last expression is valid for a sample of multĳet events, labelled by the index , and the indices , and run over all jets in a given event. Here, T = T / T is the normalised transverse energy of jet , is the angle in the transverse plane between jet and jet and ( ) is the Dirac delta function, which ensures = . The distribution is normalised to the total cross-section, .
To cancel out uncertainties that are symmetric in cos , the ATEEC function is defined as the difference between the forward (cos > 0) and the backward (cos < 0) part of the TEEC function, i.e.
1 dΣ asym d cos Both the TEEC and ATEEC functions are sensitive to gluon radiation and show a clear dependence on the strong coupling. Numerical results at next-to-leading order (NLO) for the jet-based TEEC function were obtained [43] by using the Nlojet++ program [44,45], which provides the LO and NLO calculations for three-jet production. Furthermore, numerical results for the hadron-based TEEC function at next-to-leading order plus next-to-next-to-leading-logarithmic (NLO+NNLL) accuracy were computed recently [46].
The ATLAS Collaboration has presented measurements of the TEEC and ATEEC functions at centre-ofmass energies of √ = 7 TeV [47] and √ = 8 TeV [48], determining s as a function of the interaction scale and testing the running of s predicted by the renormalisation group equation (RGE). The existence of new coloured fermions would imply modifications to the QCD -function [49,50]. Therefore, the running of s is not only important as a precision test of QCD at large scales but also as a test for new physics. An interpretation of these measurements in terms of constraints on new coloured particles through their impact on the running of s is presented in Ref. [51].
This analysis presents a measurement of the TEEC and ATEEC functions at a centre-of-mass energy of √ = 13 TeV using 139 fb −1 of collision data, recorded by the ATLAS detector. This extends the energy range of the measurement with respect to previous analyses and improves the experimental precision on these observables.
The publication of the next-to-next-to-leading-order (NNLO) corrections to three-jet production in collisions [52,53] provides a significant theoretical improvement in the precision of these observables. In particular, the theoretical uncertainties are significantly reduced with the addition of higher-order corrections. This allows more precise tests of pQCD and an important reduction of the uncertainties in the determination of the strong coupling constant .

ATLAS detector
The ATLAS detector [54] at the LHC covers nearly the entire solid angle around the collision point. 1 It consists of an inner tracking detector surrounded by a thin superconducting solenoid, electromagnetic and hadron calorimeters, and a muon spectrometer incorporating three large superconducting air-core toroidal magnets.
The inner-detector system (ID) is immersed in a 2 T axial magnetic field and provides charged-particle tracking in the range of | | < 2.5. The high-granularity silicon pixel detector covers the vertex region and typically provides four measurements per track, the first hit normally being in the insertable B-layer (IBL) installed before Run 2 [55,56]. It is followed by the silicon microstrip tracker (SCT), which usually provides eight measurements per track. These silicon detectors are complemented by the transition radiation tracker (TRT), which enables radially extended track reconstruction up to | | = 2.0. The TRT also provides electron identification information based on the fraction of hits (typically 30 in total) above a higher energy-deposit threshold corresponding to transition radiation.
The calorimeter system covers the pseudorapidity range of | | < 4.9. Within the region | | < 3.2, electromagnetic calorimetry is provided by barrel and endcap high-granularity lead/liquid-argon (LAr) calorimeters, with an additional thin LAr presampler covering | | < 1.8 to correct for energy loss in material upstream of the calorimeters. Hadron calorimetry is provided by the steel/scintillator-tile calorimeter, segmented into three barrel structures with | | < 1.7, and two copper/LAr hadron endcap calorimeters. The solid-angle coverage is completed with forward copper/LAr and tungsten/LAr calorimeter modules optimised for electromagnetic and hadronic energy measurements respectively.
The muon spectrometer (MS) comprises separate trigger and high-precision tracking chambers measuring the deflection of muons in a magnetic field generated by the superconducting air-core toroidal magnets. The field integral of the toroids ranges between 2.0 and 6.0 T m across most of the detector. Three layers of precision chambers, each consisting of layers of monitored drift tubes, covers the region | | < 2.7,

Data and Monte Carlo samples
The data used in this analysis was collected from 2015 to 2018 at a centre-of-mass energy of √ = 13 TeV. After applying quality criteria to ensure good ATLAS detector operation, the total integrated luminosity useful for data analysis is 139 fb −1 [59]. The average number of inelastic interactions produced per bunch crossing, hereafter referred to as 'pile-up', is ⟨ ⟩ = 33.6.
Several Monte Carlo (MC) simulated samples are used for this analysis; they differ in the matrix element calculation (ME) and/or the parton shower (PS). The samples were produced using the Pythia 8 [60,61], Sherpa [62,63] and Herwig 7 [64][65][66] generators. The main features of the samples are summarised in Table 1. The Pythia 8 sample was generated using Pythia 8.235. The matrix element was calculated for the 2 → 2 process. The parton shower algorithm included initial-and final-state radiation based on the dipole-style T -ordered evolution, including →¯branchings and a detailed treatment of the colour connections between partons [60]. The renormalisation and factorisation scales were set to the geometric mean of the squared transverse masses of the two outgoing particles (labelled 3 and 4), i.e.
. The NNPDF 2.3 LO parton distribution function (PDF) set [67] was used in the ME generation, in the parton shower, and in the simulation of multi-parton interactions (MPI). The ATLAS A14 [68] set of tuned parameters (tune) was used for the parton shower and MPI, while hadronisation was modelled using the Lund string model [69,70].
The Sherpa sample was generated using Sherpa 2.1.1. The matrix element calculation includes 2 → 2 and 2 → 3 processes at LO, merged using the CKKW method [71]. The Sherpa parton shower [72,73] with T ordering was used for the quark and gluon emissions. The matrix element renormalisation and factorisation scales for 2 → 2 processes were set to the harmonic mean of the Mandelstam variables , and [74], whereas the Catani-Marchesini-Webber (CMW) [75] scale was chosen for 2 → 3 processes. The CT14 NNLO [76] PDF set was used for the matrix element calculation, while the parameters used for the modelling of the MPI and the parton shower were set according to the CT10 tune [77]. The Sherpa sample made use of the dedicated Sherpa AHADIC model for hadronisation [78], which is based on the cluster fragmentation algorithm [79].
Finally, two Herwig 7 samples were generated using Herwig 7.1.3 at next-to-leading order for the 2 → 2 process. The matrix element was calculated using Matchbox [80] with the MMHT2014 NLO PDF [81]. The renormalisation and factorisation scales were set to the T of the leading jet. The first sample used an angle-ordered parton shower, while the second sample used a dipole-based parton shower. In both cases, the parton shower was interfaced to the matrix element calculation using the MC@NLO matching scheme [82]. The angle-ordered shower evolves on the basis of 1→2 splittings using a generalised angular variable and employs a global recoil scheme once showering has terminated. The dipole-based shower uses 2→3 splittings with Catani-Seymour kernels with an ordering in transverse momentum that enables it to perform recoils on an emission-by-emission basis. For both Herwig 7 samples, the parameters that control the MPI and parton shower simulation were set according to the MMHT2014 tune [81], and the hadronisation was modelled by the cluster fragmentation algorithm.
The Pythia 8, Sherpa and Herwig 7 samples were passed through the Geant4-based [83] ATLAS detector-simulation programs [84], since they were used to unfold the measurements to the particle level as described in Section 5. They were reconstructed and analysed with the same processing chain as the data. The generation of the simulated event samples includes the effect of multiple interactions per bunch crossing, simulated using the Pythia A3 tune [85], and the effect of interactions from bunch crossings before or after the one containing the hard interaction.
| | < 2.4 to reject pile-up jets and reduce experimental uncertainties. In addition, jets are required to satisfy quality criteria that reject non-collision background and noise jets [91]. More than 99% of all signal jets pass these criteria. Events are required to have at least two selected jets. The two leading jets are further required to have the scalar sum of their transverse momenta, T2 = T1 + T2 > 1 TeV. This ensures a trigger efficiency of ≈100%. About 57.5 million events in data satisfy the selection criteria. The data are binned in ten intervals of T2 in which the TEEC and ATEEC distributions are measured. The MC predictions are reweighted using the distribution of the average number of interactions per bunch crossing, so that the pile-up conditions in the simulation exactly match those in data.

Unfolding to particle level
To make detector-independent comparisons with particle-level parton shower MC predictions and fixedorder theoretical predictions, the measured distributions are corrected for distortions induced by the response of the ATLAS detector and associated reconstruction algorithms. The fiducial phase-space region is defined at particle level for all particles with a decay length > 10 mm; these particles are referred to as 'stable'. Particle-level jets are reconstructed using the anti-algorithm with = 0.4 using stable particles, excluding muons and neutrinos. The fiducial phase space closely follows the event selection criteria defined in Section 4. Particle-level jets are required to have T > 60 GeV and | | < 2.4. Events with at least two particle-level jets are considered. These events are also required to have T2 > 1 TeV.
The unfolding is performed using an iterative algorithm based on Bayes theorem [92], which takes into account inefficiencies, acceptance and resolution effects in the detector response that lead to bin migrations between the detector-level and particle-level phase spaces. For each observable, the method makes use of a transfer matrix, , obtained from MC simulation, that parameterises the probability of a jet pair generated in bin being reconstructed in bin . The binning of the distributions is chosen as a compromise between maximising the number of bins while minimising migration between bins due to the resolution of the measured variables. This ensures that the purity from migrations is above 60% in all regions of the phase space. The inverse problem can be written as a linear equation The quantities and are the contents of the detector-level distribution in bin and the contents of the particle-level distribution in bin , respectively, while the factors E and P are efficiency and purity corrections, which are estimated by using the MC simulation.
The unfolding is performed separately in each bin of T2 only considering the dependence on cos . The effect of migrations across different T2 bins is small and is treated by the efficiency and purity corrections. It has been verified that using a two-dimensional transfer matrix, including T2 and cos , gives results compatible within the statistical uncertainties.
Five iterations are necessary for the procedure to converge, i.e. for the relative difference between one iteration and the next, averaged over all bins of the distribution, to be below one per mille. The statistical uncertainties in both data and MC samples are propagated through the unfolding procedure by using pseudo-experiments. A set of 1000 replicas is considered for each measured distribution by applying a Poisson-distributed fluctuation around the nominal measured distribution. Each of these replicas is unfolded using a fluctuated version of the transfer matrix, which produces the corresponding set of 1000 replicas of the unfolded spectra. The statistical uncertainty is defined as the standard deviation of all replicas.
The results obtained by unfolding the data with Pythia 8 are used to obtain the nominal differential crosssections, whereas Sherpa and Herwig 7 MC predictions are used to estimate the systematic uncertainty due to the model dependence, as discussed in Section 6.

Experimental uncertainties
The dominant sources of systematic uncertainty in the measurements arise from knowledge of the jet energy scale and resolution and the modelling of the final state. The detector-level distributions are recomputed and unfolded for each systematic variation. The systematic uncertainties are then propagated through the unfolding via their impact on the transfer matrices, the purity and the efficiency. A detailed description of the systematic uncertainties and their typical values for the TEEC and ATEEC functions is given below: • Jet Energy Scale and Resolution: The jet energy scale (JES) and jet energy resolution (JER) uncertainties are estimated as described in Ref. [89,90]. The JES uncertainties consist of a set of 113 independent components, which depend on the jet T and . These components account for different effects, including the different response to quarks and gluons and the flavour composition of the sample. Sources related to the in situ methods for the determination of JES [89, 90] are also included. The total JES uncertainty in the T value of individual jets varies from approximately 2% at T = 60 GeV to below 1% for T > 500 GeV, with a mild dependence on . The JER uncertainty is estimated by using 34 independent components, arising from the determination of the JER noise term and the different in situ methods [90]. The effect of the total JER uncertainty is evaluated by smearing the energy of the jets in the MC simulation by about 2% at T = 60 GeV to about 0.5% for T of several hundred GeV. Each JES and JER uncertainty component is independently propagated by varying the energy and T of each jet by one standard deviation. The total uncertainty in the TEEC distribution ranges from a maximum of 2% for 1 TeV < T2 < 1.2 TeV to 1.5% for T2 > 3.5 TeV while, for the ATEEC, it is below 1% in the whole phase space of the analysis.
• Jet Angular Resolution: The jet angular resolution (JAR) uncertainty is conservatively estimated by smearing the azimuthal coordinate of the jets by their resolution in MC simulation. The variations are performed with the T component of the jets held constant. The value of the JAR uncertainty is below 0.5% for both the TEEC and ATEEC distributions in all regions of phase space.
• Unfolding: The mismodelling of the data by the MC simulation is accounted for as an additional source of uncertainty. This is assessed by reweighting the particle-level distributions so that the detector-level distributions predicted by the MC samples match those in the data. The modified detector-level MC distributions are then unfolded using the method described in Section 5. The difference between the modified particle-level MC distribution and the original particle-level MC distribution is taken as the uncertainty. This uncertainty is negligible for both the TEEC and ATEEC in each of the phase space regions of the analysis.
• Monte Carlo Modelling: The modelling uncertainty is estimated by comparing the distributions obtained by using Pythia 8, Sherpa and Herwig 7 in the unfolding procedure. The envelope of the differences between the estimated cross-sections defines the systematic uncertainty. The value of this uncertainty for the TEEC varies from 1.5% to 2%, depending on the phase space region while, for the ATEEC, it is always below 0.5% and has an impact only for cos < −0.95.
A breakdown of the relative systematic uncertainties in the measurement of the TEEC and ATEEC distributions is shown in Figure 1.  The total systematic uncertainty is estimated by adding the effects previously listed in quadrature. In addition, the statistical uncertainty in the data and MC simulation is propagated to the differential cross-sections through the unfolding procedure using pseudo-experiments to properly take into account the statistical correlations between bins. The pseudo-experiments are also used to estimate the statistical component of each systematic uncertainty. This statistical component is mitigated using the Gaussian Kernel smoothing technique [93]. The values provided above are quoted after application of this procedure. The impact of these sources of uncertainty depends on the region of phase space. Overall, the JES uncertainty dominates both the TEEC and ATEEC measurements, while the MC modelling and JER uncertainties are also important for the TEEC and ATEEC measurements, respectively.

Experimental results
The unfolded TEEC and ATEEC distributions are presented and compared with the MC predictions described in Section 3. Figure 2 shows this comparison for the TEEC and ATEEC distributions obtained in the inclusive region T2 > 1 TeV. Additionally, Figure 3 shows this comparison for both the TEEC and ATEEC for two different bins of T2 .   Figure 2: Comparison between the data and the MC predictions for (a) the TEEC and (b) the ATEEC for the inclusive T2 bin. The data error bands contain both statistical and systematic components, summed in quadrature. The lower panels show the ratios of the MC predictions to the data, including the statistical uncertainty on the MC predictions. The value of the first bin of the ATEEC distribution is negative and, therefore, not represented in logarithmic scale.     The TEEC and ATEEC distributions presented in Figures 2 and 3 show the features observed in previous results [47,48]. The TEEC shows two peaks arising from back-to-back (cos = −1) and collinear (cos = 1) configurations, together with a central plateau arising from the wide-angle radiation. The effect of the jet radius is seen as a kink in the TEEC distributions at cos ≃ 0.92. The ATEEC exhibits a steep fall-off of several orders of magnitude between cos = −1 and cos = 0.
At low values of T2 , the best description of the TEEC distributions is achieved by both Sherpa and Herwig 7 with the angle-ordered parton shower. The dipole parton shower in Herwig 7 does not describe the data well, in particular in the region | cos | > 0.4, where discrepancies of up to 20% are observed close to the edges of the TEEC distribution. Pythia 8 shows discrepancies compared with data for −0.7 < cos < −0.2, where differences of O (5% − 10%) are observed. A good description of the ATEEC distribution is achieved by both Sherpa and Herwig 7, while Pythia 8 underestimates the measured distribution for cos > −0.7 and the Herwig 7 sample matched to the dipole parton shower underestimates the values of the ATEEC for cos < −0.7.
For large values of T2 , Pythia 8 gives the best description of the TEEC data, while both Sherpa and Herwig 7 tend to overestimate the height of the central plateau. The dipole parton shower in Herwig 7 shows a similar behaviour to the low T2 region, with the larger discrepancies observed near the edges of the TEEC distribution, although limited to the region | cos | > 0.7. For the ATEEC, Pythia 8 and Sherpa give a good description of the data, while the Herwig 7 sample matched to the angle-ordered parton shower overestimates the values of the ATEEC for cos < −0.7. The cancellation of discrepancies in -symmetric regions between the data and the dipole shower in Herwig 7 brings a reasonable agreement for the ATEEC in the high energy regime.

Theoretical predictions and uncertainties
The theoretical predictions for the TEEC and ATEEC functions are calculated using pQCD at LO, NLO and NNLO in powers of the strong coupling constant, s ( R ), as discussed in Ref. [52,53]. The calculations make use of OpenLoops2 [94], FivePointAmplitudes [95] and PentagonFunctions++ [96] for the evaluation of the scattering amplitudes. A mixed flavour scheme is assumed. On the one hand, the effect of the top quark is taken into account when evaluating the radiative corrections to the gluon propagator. Thus, the solution to the RGE includes six active flavours by considering different values of Λ QCD when crossing the corresponding mass thresholds. On the other hand, the perturbative corrections to the three-jet cross-section assume f = 5 massless quark flavours, thus neglecting the effect of the top quark in both the virtual and real contributions. The contribution from real top-quark production is estimated to be below 0.3% for the phase space relevant to this analysis.
At O ( 5 ) the numerator in Eq. 1 entails the calculation of 2 → 3 subprocesses at NNLO, 2 → 4 at NLO and 2 → 5 at LO. In order to avoid the triple-collinear singularities appearing in the latter, the angular range is restricted to | cos | < 0.92. This avoids calculating the three-loop virtual corrections to the 2 → 2 subprocesses. Thus, with this azimuthal cut, the denominator in Eq. 1 includes the 2 → 2, 2 → 3 and 2 → 4 subprocesses up to, and including, O ( 4 ) corrections.
The value of s ( R ) that enters in the partonic matrix-element calculation is obtained by evolving the same s ( ) value as used to determine the PDF set under consideration. The NNLO PDF sets MMHT2014 [81], NNPDF 3.0 [97] and CT14 [76], as provided in the LHAPDF package [98], are convolved with the partonic cross-sections.
The NNLO corrections are calculated using O (10 13 ) events. They include real-real, real-virtual and virtual-virtual finite terms, as well as single-and double-unresolved terms and the finite remainder [99,100]. For illustrative purposes, Figure 4 shows sample diagrams corresponding to real-real, real-virtual and virtual-virtual contributions to the → jjj process at NNLO. The pQCD predictions for the TEEC and ATEEC functions are obtained using jets clustered from final-state partons as inputs to the anti-algorithm with = 0.4 [86]. Therefore, non-perturbative correction factors are derived to account for hadronisation and underlying event effects on the observables. These correction factors are calculated as the ratio of the MC predictions including both hadronisation and underlying event effects to the MC predictions in which these effects are not included. The central values for these correction factors are calculated using the A14 tune in Pythia 8 [68], and are close to unity to within 0.5%, except for cos > 0.86 where they go down to as far as 0.95.
The uncertainties in the theoretical predictions arise from three sources: contributions from higher-order corrections estimated using scale variations, the uncertainties in the PDF and the uncertainties in the non-perturbative corrections to the pQCD predictions.
• The scale uncertainties are evaluated by considering independent variations of the renormalisation and factorisation scales by a factor of two, excluding those in which R and F are varied in opposite directions. As in previous analyses [47,48], the scale variations are treated as correlated between the three-jet and two-jet cross sections, respectively in the numerator and denominator of Eq. 1. The uncertainty in the distributions is defined as the envelope of the six resulting variations. Overall, the scale uncertainty bands are asymmetric and larger for the lower part of the TEEC and ATEEC distributions. The value of the scale uncertainties are below 2% in all regions of phase space, while they account for up to 6% at NLO. Thus, this source of uncertainty is reduced by a factor of three by the NNLO corrections. This is the dominant theoretical uncertainty for both TEEC and ATEEC.
• The PDF uncertainties are obtained by considering the set of eigenvectors / replicas provided by each PDF group [76,81,97], and are propagated to the predicted distributions using the corresponding recommendations from each particular PDF set. Due to the large computing time required for the NNLO calculations, the PDF variations are calculated at NLO and extrapolated to the NNLO predictions using their relative values. The resulting uncertainties are below 1% for the low T2 bins, and can reach up to 2% for the higher T2 bins.
• The uncertainties in the non-perturbative corrections are evaluated by considering the envelope of the differences between five different predictions and the Pythia 8 A14 tune [68]. These predictions were generated using Pythia 8.235 and Herwig 7.2.1 with the angle-ordered parton shower. They include the 4C [102] and Monash [103] Pythia 8 tunes, and the Herwig 7 default, soft and baryonic reconnection tunes [104]. The value of this uncertainty is below 1% in the phase-space region | cos | < 0.92, where s is determined.
The relative values of the theoretical uncertainties for the TEEC and ATEEC distributions are shown in Figure 5 for the inclusive T2 bin. The effect of varying ( ) within the range 0.117 -0.119 is also shown for illustrative purposes. (a)  Figures 6 and 7 show the ratios of the data, as well as of the LO and NNLO theoretical predictions to the NLO predictions for the TEEC and ATEEC, respectively. The description of the data provided by the NNLO pQCD predictions is very good, improving with respect to the NLO prediction. In these figures, the distributions are rebinned to reduce statistical fluctuations in the theoretical predictions that, despite the high statistics, can be present due to insufficient precision on the cancellation of the different terms. The reduction of the scale uncertainties is made evident from these figures, as well as the improvement in the description. In particular, the region with cos > 0.80, shows a large improvement for all T2 bins. For the higher T2 bins, and for the MMHT2014 PDF, the NNLO theoretical predictions are slightly above the data. Since the CT14 PDF shows a better behaviour in this region of the phase space, this effect may be due to the limited accuracy of the PDF determinations at high values of the Bjorken variable . The results obtained using the CT18 PDF set agree well with those for MMHT2014.

Determination of the strong coupling constant
The strong coupling constant at the scale of the pole mass of the boson, s ( ) can be determined from the comparison of the data with the theoretical predictions by considering the following 2 function where the theoretical predictions are varied according to In Eq.
(2) and Eq. (3), s stands for s ( ); is the value of the -th point of the distribution as measured in data, while Δ is its statistical uncertainty. The statistical uncertainty in the theoretical predictions dominates over the statistical uncertainty in data, and is also included as Δ , while ( ) is the relative value of the -th correlated source of systematic uncertainty in bin .
This technique takes into account the correlations between the different sources of systematic uncertainty discussed in Section 6 by introducing the nuisance parameters { }, one for each source of experimental uncertainty, each of which follow a standard normal distribution. Thus, the minimum of the 2 function defined in Eq. (2) is found in a 150-dimensional space, in which 149 correspond to nuisance parameters { } and one to s ( ). The uncertainty due to the choice of the MC model in the unfolding is treated by performing alternative fits in which the data distribution is unfolded using different MC models, including Herwig 7 and Sherpa. The envelope of the values obtained for each model is used as the systematic uncertainty.
The method also requires an analytical expression for the dependence of each observable on the strong coupling constant, which is given by ( s ) for bin . For each PDF set, the corresponding s ( ) variation range is considered and the theoretical prediction is obtained for each value of s ( ). The functions ( s ) are then obtained by fitting the predicted values of the TEEC (ATEEC) in each ( T2 , cos ) bin to a third-order polynomial in .
For both the TEEC and ATEEC functions, the fits to extract s ( ) are repeated separately for each T2 interval, thus determining a value of s ( ) for each energy bin. The theoretical uncertainties are determined by shifting the theory distributions by each of the uncertainties separately, recalculating the functions ( s ) and determining a new value of s ( ). The uncertainty is determined by taking the difference between this value and the nominal one. The uncertainty due to the choice of the MC tune in the non-perturbative correction is estimated by repeating the fit with the theory distributions corrected with the other tunes, and taking the envelope relative to the nominal fit.
Each of the fitted values of s ( ), which are obtained in the MS subtraction scheme, is then evolved to the corresponding measured scale using the NNLO solution to the RGE. When evolving s ( ) to s ( ), the appropriate matching conditions for the strong coupling constant at the f = 5 and f = 6 flavour thresholds are applied so that s ( ) is a continuous function across quark thresholds. The value of is obtained as the average value of =ˆT calculated at NNLO for each T2 bin.
The values of ( ) obtained from global fits to the TEEC and ATEEC distributions, using the MMHT 2014, CT14 and NNPDF 3.0 PDF sets are shown in Tables 2 and 3, together with the corresponding 2 values. The values of s ( ) determined from fits to the TEEC functions, both in the inclusive and exclusive T2 bins, are presented in Table 4, while the values obtained from fits to the ATEEC functions are presented in Table 5. The average values and the reduced 2 values at the minimum are also shown. The values displayed in Tables 4 and 5 are obtained using the MMHT 2014 PDF in the theoretical predictions. Table 2: Values of ( ) obtained from fits to the TEEC distributions to the NNLO theoretical predictions obtained using different PDF sets. The uncertainty referred to as represents the theoretical scale uncertainty. The uncertainty referred to as NP is related to the non-pQCD corrections. The uncertainty referred to as (mod.) is related to the MC model used in the unfolding. The value of the 2 function at its minimum, as well as the number of degrees of freedom is also indicated.  Table 3: Values of ( ) obtained from fits to the ATEEC distributions to the NNLO theoretical predictions obtained using different PDF sets. The uncertainty referred to as represents the theoretical scale uncertainty. The uncertainty referred to as NP is related to the non-pQCD corrections. The uncertainty referred to as (mod.) is related to the MC model used in the unfolding. The value of the 2 function at its minimum, as well as the number of degrees of freedom is also indicated. The central values of , obtained with different PDFs, in Tables 2 and 3 show differences that are, in some cases, larger than the uncertainty band on the eigenvectors or replicas provided for each PDF set. This is understood to arise from the differences between PDFs at high values of the Bjorken variable , which influence the theory description of the TEEC and ATEEC in the phase space region relevant for this analysis.
Since The central values of determined from both observables are correlated with a Pearson correlation coefficient = 0.86 ± 0.02 (exp.), and are found to be compatible within the quoted uncertainties. While the value determined from the fit to the TEEC distributions has a better experimental precision, mainly due to its smaller statistical uncertainty, the value determined from the fit to the ATEEC distributions exhibits a better theoretical precision.
The fitted values of s ( ) for each bin, together with the values of the nuisance parameters, are used to evaluate the agreement between the data and the fitted predictions. Figure 8 shows the ratio of the TEEC data to the fitted theoretical distributions for both the inclusive and exclusive bins, together with the theoretical and experimental uncertainties. Similarly, Figure 9 shows the ratios of data to theory for the ATEEC distributions. An excellent agreement is observed for both the TEEC and ATEEC measurements in all regions of phase space.
The values of s ( ) obtained in Tables 4 and 5 are evolved to the corresponding energy scale using the RGE as explained above. This allows testing the asymptotic behaviour of QCD by comparing the measured points with the prediction from the RGE. Figures 10 and 11 show the values of s ( ) together with the world average band provided by the Particle Data Group [105] and values of s obtained in other analyses [47,48,[106][107][108][109][110][111][112]. To compare the results with other experiments, the value of is chosen as half of the average value ofˆT for each T2 bin. The results show a good agreement between all measurements and the RGE prediction. The values of obtained from fits to the theoretical predictions using different PDF sets, shown in Figs. 12 and 13, are found to be in good agreement. Since the main purpose of the present analysis is to provide a first determination of using NNLO predictions for three-jet cross sections, a combined determination of and the PDF [113] is not attempted. Such analysis would require to use very large computing resources and is well out of the scope of this paper.

Conclusions
A measurement of transverse energy-energy correlations and their corresponding asymmetries in multĳet events using 139 fb −1 of collisions at √ = 13 TeV collected with the ATLAS detector at the LHC is presented. High-energy multĳet events are selected by requiring the scalar sum of T of the two leading jets, T2 , to be above 1 TeV, and the data are binned in this variable to study the scale dependence of these observables.
The data are corrected for detector effects and systematic uncertainties are evaluated. The experimental uncertainties are dominated by the uncertainty in the jet energy scale and the model used to correct for detector effects. The total uncertainty is found to be of the order of 2% for the TEEC and 1% for the ATEEC. The results are compared with MC predictions from different generators, including Pythia 8, Sherpa and Herwig 7 with two different parton showers, one angle-ordered and one based on dipole radiation.
The data are also compared with novel theoretical predictions at next-to-next-to-leading order in perturbative QCD, which are corrected for non-perturbative effects such as hadronisation and multiple parton interactions. The theoretical uncertainty due to the choice of the renormalisation and factorisation scales amount to 2% for the TEEC and 3% for the ATEEC, being the dominant uncertainty for both observables. The agreement between the data and the theoretical predictions is excellent. for the TEEC and ATEEC, respectively. The inclusion of next-to-next-to-leading-order corrections to three-jet production in the theoretical predictions for the first time has allowed a reduction by a factor of 3 of the theoretical uncertainties in both the cross-section calculation for the TEEC and ATEEC distributions and in the determination of the strong coupling constant .    [114] ATLAS Collaboration, ATLAS Computing Acknowledgements, ATL-SOFT-PUB-2021-003, 2021, url: https://cds.cern.ch/record/2776662.