\section{Introduction}
Multijet final states, produced in proton--proton ($pp$) collisions with large momentum transfer at the LHC provide an ideal testing ground for perturbative Quantum Chromodynamics (pQCD). Event shapes~\cite{evShapes1, evShapes2} 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, $\alpha_{\mathrm{s}}$. Event shape variables were measured in $e^+ e^-$ collision experiments from PETRA and PEP~\cite{plutoEventShapes,tassoEventShapes,jadeEventShapes} to LEP and SLC~\cite{opalEventShapes,alephEventShapes,thrustDelphi,l3eventShapes,SLD}, at the $ep$ collider HERA~\cite{zeusEventShapes1,zeusEventShapes2,h1EventShapes1,h1EventShapes2,h1EventShapes3} and in hadron--hadron collisions at the Tevatron~\cite{Aaltonen} and the LHC~\cite{Khachatryan, Aad, evShapes}.
 
A particularly interesting, infrared safe, event-shape observable is the energy--energy correlation (EEC) function, which was originally introduced for $e^{+} e^{-}$ annihilation experiments~\cite{EEC-TH-00, EEC-TH-01}. The EEC function and its associated azimuthal asymmetry (AEEC) can be calculated in pQCD, and the $\mathcal{O}(\alpha^{2}_{\mathrm{s}})$ corrections are modest~\cite{EEC-TH-02, EEC-TH-03, EEC-TH-04, EEC-TH-04b, DESY-88-081}. They were also studied in the nearly back-to-back limit~\cite{EEC-TH-05}. The EEC measurements~\cite{PLUTO, MARKII, MARKJ, CELLO, JADE, MAC, TASSO, TOPAZ, DELPHI, OPAL, ALEPH, L3, eecSLD} have had significant impact on the early precision tests of QCD and in the determination of the strong coupling constant.
 
Transverse energy--energy correlations (TEEC) and their associated azimuthal asymmetries (ATEEC) were proposed as the appropriate generalisation for hadron collider experiments in Ref.~\cite{CERN-TH-3800}, 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 $E_{\mathrm{T}} = E / \cosh y$, where $E$ is the jet energy and $y$ 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~\cite{TEEC-TH-00}, i.e.
\[
\frac{1}{\sigma}\frac{\mathrm{d}\Sigma}{\mathrm{d}\cos\phi} \equiv \frac{1}{\sigma}\sum_{ij}\int\frac{\mathrm{d}\sigma}{\mathrm{d}x_{\mathrm{T}i} \mathrm{d}x_{\mathrm{T}j} \mathrm{d}\cos\phi}x_{\mathrm{T}i}x_{\mathrm{T}j} \mathrm{d}x_{\mathrm{T}i} \mathrm{d}x_{\mathrm{T}j} = \frac{1}{N} \sum_{A=1}^N \sum_{ij} \frac{E^A_{\mathrm{T}i}E^A_{\mathrm{T}j}}{\left(\sum_k E^A_{\mathrm{T}k}\right)^2}\delta(\cos \phi - \cos\varphi_{ij}),
\]
where the last expression is valid for a sample of $N$ multijet events, labelled by the index $A$, and the indices $i$, $j$ and $k$ run over all jets in a given event. Here, $x_{\mathrm{T}i} = E_{\mathrm{T}i} / \sum_k E_{\mathrm{T}k}$ is the normalised transverse energy of jet $i$, $\varphi_{ij}$ is the angle in the transverse plane between jet $i$ and jet $j$ and $\delta(x)$ is the Dirac delta function, which ensures $\phi = \varphi_{ij}$. The normalisation to the total dijet cross-section, $\sigma$, ensures that the integral of the TEEC function over $\cos\phi$ is unity.
 
To cancel out uncertainties that are symmetric in $\cos\phi$, the ATEEC function is defined as the difference between the forward ($\cos\phi > 0$) and the backward ($\cos\phi < 0$) part of the TEEC function, i.e.
\[
\frac{1}{\sigma} \frac{\mathrm{d}\Sigma^{\mathrm{asym}}}{\mathrm{d}\cos\phi}  =
\frac{1}{\sigma} \frac{\mathrm{d}\Sigma}{\mathrm{d}\cos\phi} \biggr\rvert_\phi -
\frac{1}{\sigma} \frac{\mathrm{d}\Sigma}{\mathrm{d}\cos\phi} \biggr\rvert_{\pi-\phi} .
\]
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~\cite{TEEC-TH-01} by using the \textsc{Nlojet++} program~\cite{IPPP/01/48, TEEC-TH-03}, 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~\cite{TEEC-TH-02}.
 
The ATLAS Collaboration has presented measurements of the TEEC and ATEEC functions at centre-of-mass energies of $\sqrt{s} = 7~\TeV$~\cite{STDM-2014-03} and $\sqrt{s} = 8~\TeV$~\cite{STDM-2016-10}, determining $\alpha_{\mathrm{s}}$ as a function of the interaction scale and testing the running of $\alpha_{\mathrm{s}}$ predicted by the renormalisation group equation (RGE). The existence of new coloured fermions would imply modifications to the QCD $\beta$-function~\cite{TEEC-TH-04, TEEC-TH-05}. Therefore, the running of $\alpha_{\mathrm{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 $\alpha_{\mathrm{s}}$ is presented in Ref.~\cite{TEEC-TH-06}.
 
This analysis presents a measurement of the TEEC and ATEEC functions at a centre-of-mass energy of $\sqrt{s} = 13~\TeV$ using 139~\ifb of $pp$ 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 $pp$ collisions~\cite{mitov2022,mitov2023} 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 $\alpha_s$.


\section{ATLAS detector}
\newcommand{ATLAS uses a right-handed coordinate system with its origin at the nominal interaction point (IP) in the centre of the detector and the \(z\)-axis along the beam pipe. The \(x\)-axis points from the IP to the centre of the LHC ring, and the \(y\)-axis points upwards. Cylindrical coordinates \((r,\phi)\) are used in the transverse plane, \(\phi\) being the azimuthal angle around the \(z\)-axis. The pseudorapidity is defined in terms of the polar angle \(\theta\) as \(\eta = -\ln \tan(\theta/2)\).}{ATLAS uses a right-handed coordinate system with its origin at the nominal interaction point (IP) in the centre of the detector and the \(z\)-axis along the beam pipe. The \(x\)-axis points from the IP to the centre of the LHC ring, and the \(y\)-axis points upwards. Cylindrical coordinates \((r,\phi)\) are used in the transverse plane, \(\phi\) being the azimuthal angle around the \(z\)-axis. The pseudorapidity is defined in terms of the polar angle \(\theta\) as \(\eta = -\ln \tan(\theta/2)\).}
 
The ATLAS detector~\cite{PERF-2007-01} at the LHC covers nearly the entire solid angle around the collision point.\footnote{ATLAS uses a right-handed coordinate system with its origin at the nominal interaction point (IP) in the centre of the detector and the \(z\)-axis along the beam pipe. The \(x\)-axis points from the IP to the centre of the LHC ring, and the \(y\)-axis points upwards. Cylindrical coordinates \((r,\phi)\) are used in the transverse plane, \(\phi\) being the azimuthal angle around the \(z\)-axis. The pseudorapidity is defined in terms of the polar angle \(\theta\) as \(\eta = -\ln \tan(\theta/2)\).}
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 \SI{2}{\tesla} axial magnetic field
and provides charged-particle tracking in the range of \(|\eta| < 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~\cite{ATLAS-TDR-19,PIX-2018-001}.
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 \(|\eta| = 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 \(|\eta| < 4.9\).
Within the region \(|\eta|< 3.2\), electromagnetic calorimetry is provided by barrel and
endcap high-granularity lead/liquid-argon (LAr) calorimeters,
with an additional thin LAr presampler covering \(|\eta| < 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 \(|\eta| < 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 \num{2.0} and \SI{6.0}{\tesla\metre}
across most of the detector.
Three layers of precision chambers, each consisting of layers of monitored drift tubes, covers the region \(|\eta| < 2.7\),
complemented by cathode-strip chambers in the forward region, where the background is highest.
The muon trigger system covers the range of \(|\eta| < 2.4\) with resistive-plate chambers in the barrel, and thin-gap chambers in the endcap regions.
 
Interesting events are selected by the first-level trigger system implemented in custom hardware,
followed by selections made by algorithms implemented in software in the high-level trigger~\cite{TRIG-2016-01}.
The first-level trigger accepts events from the \SI{40}{\MHz} bunch crossings at a rate below \SI{100}{\kHz},
which the high-level trigger further reduces to record events to disk at about \SI{1}{\kHz}. An extensive software suite~\cite{ATL-SOFT-PUB-2021-001} is used in data simulation, in the reconstruction and analysis of real and simulated data, in detector operations, and in the trigger and data acquisition
systems of the experiment.


\section{Data and Monte Carlo samples}
\label{sec:MC_simulation}
The data used in this analysis was collected from 2015 to 2018 at a centre-of-mass energy of $\sqrt{s} = 13~\TeV$. After applying quality criteria to ensure good ATLAS detector operation, the total integrated luminosity useful for data analysis is 139~$\ifb$~\cite{DAPR-2018-01}. The average number of inelastic $pp$ interactions produced per bunch crossing, hereafter referred to as `pile-up', is $\langle\mu\rangle = 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 \textsc{Pythia 8} \cite{Sjostrand:2006za, Sjostrand:2014zea}, \textsc{Sherpa} \cite{sherpa,sherpa22} and \textsc{Herwig 7} \cite{Bahr:2008pv,Bellm:2015jjp,Bellm:2017bvx} generators. The main features of the samples are summarised in Table~\ref{tab:mcSamples}.
\begin{table}[H]
\small
\caption{Properties of the MC samples used in the analysis, including the perturbative order in $\alpha_{\mathrm{s}}$, the number of final-state partons, the PDF set, the parton shower algorithm, the renormalisation and factorisation scales $\mu_R$ and $\mu_F$, and the value of $\alpha_{\mathrm{s}}(m_\Zboson)$ for the matrix element. \textsc{Herwig 7} was used for the generation of two samples, differing by their treatment of the parton shower.}
\label{tab:mcSamples}
\begin{center}
\begin{tabular}{c|c|c|c|c|c|c}
\T\B Generator & ME order & ME partons & PDF set & Parton shower & Scales $\mu_R,\mu_F$ & $\alpha_{\mathrm{s}}(m_\Zboson)$\\
\hline
\hline
\T\T\B \textsc{Pythia 8} & LO & 2 & \textsc{NNPDF 2.3 LO} & $\pt$-ordered & $(m^2_{\mathrm{T}3}\cdot m^2_{\mathrm{T}4})^\frac{1}{2}$ & 0.140\\
\hline
\multirow{2}{*}{\T\B \textsc{Sherpa}} & \multirow{2}{*}{LO} & \multirow{2}{*}{2,3} & \multirow{2}{*}{\textsc{CT14 NNLO}} & \multirow{2}{*}{CSS (dipole)} & \T\B $H(s,t,u)$ [$2\to 2$] & \multirow{2}{*}{0.118}\\
&    &     &                    &              & CMW [$2\to 3$]                       & \\
\hline
\multirow{2}{*}{\T\B \textsc{Herwig 7}} & \multirow{2}{*}{NLO} & \multirow{2}{*}{2,3} & \multirow{2}{*}{\textsc{MMHT2014 NLO}} & \T\B Angular-ordered & \multirow{2}{*}{$\max\left\{p_{\mathrm{T}i}\right\}_{i=1}^N$} & \multirow{2}{*}{0.120}\\
&                                    &  & & Dipole &\\
\end{tabular}
\end{center}
\end{table}
The \textsc{Pythia 8} sample was generated using \textsc{Pythia 8.235}. The matrix element was calculated for the $2\to 2$ process. The parton shower algorithm included initial- and final-state radiation based on the dipole-style $\pt$-ordered evolution, including $\gamma\to q\bar{q}$ branchings and a detailed treatment of the colour connections between partons~\cite{Sjostrand:2006za}. 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. $\sqrt{m^2_{\mathrm{T}3} \cdot m^2_{\mathrm{T}4}}$ = $\sqrt{(\pt^2 + m_3^2) \cdot (\pt^2 + m_4^2)}$. The \textsc{NNPDF 2.3 LO} parton distribution function (PDF) set~\cite{nnpdf} was used in the ME generation, in the parton shower, and in the simulation of multi-parton interactions (MPI). The ATLAS A14~\cite{ATL-PHYS-PUB-2014-021} set of tuned parameters (tune) was used for the parton shower and MPI, while hadronisation was modelled using the Lund string model~\cite{Andersson:1983ia,Sjostrand:1984ic}.
 
The \textsc{Sherpa} sample was generated using \textsc{Sherpa 2.1.1}. The matrix element calculation includes $2\to 2$ and $2\to 3$ processes at LO, merged using the CKKW method~\cite{Hoeche:2012yf}. The \textsc{Sherpa} parton shower~\cite{Catani:2001cc, Krauss:2002up} with $\pt$ ordering was used for the quark and gluon emissions. The matrix element renormalisation and factorisation scales for $2\to 2$ processes were set to the harmonic mean of the Mandelstam variables $s$, $t$ and $u$~\cite{mandelstam}, whereas the Catani--Marchesini--Webber (CMW)~\cite{CMWscale} scale was chosen for $2\to 3$ processes. The \textsc{CT14 NNLO}~\cite{ct14} 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 \textsc{CT10} tune~\cite{ct10}. The \textsc{Sherpa} sample made use of the dedicated Sherpa AHADIC model for hadronisation~\cite{Winter:2003tt}, which is based on the cluster fragmentation algorithm~\cite{Webber:1983if}.
 
Finally, two \textsc{Herwig 7} samples were generated using \textsc{Herwig 7.1.3} at next-to-leading order for the $2\to 2$ process. The matrix element was calculated using \textsc{Matchbox} \cite{Platzer:2011bc} with the \textsc{MMHT2014 NLO} PDF~\cite{mmht}. The renormalisation and factorisation scales were set to the $\pt$ 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 \textsc{MC@NLO} matching scheme~\cite{Frixione:2002ik}. The angle-ordered shower evolves on the basis of 1$\to $2  splittings using  a  generalised  angular  variable  and  employs  a  global recoil scheme once showering has terminated. The dipole-based shower uses 2$\to $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 \textsc{Herwig 7} samples, the parameters that control the MPI and parton shower simulation were set according to the \textsc{MMHT2014} tune~\cite{mmht}, and the hadronisation was modelled by the cluster fragmentation algorithm.
 
The \textsc{Pythia 8}, \textsc{Sherpa} and \textsc{Herwig 7} samples were passed through the \textsc{Geant4}-based~\cite{Agostinelli:2002hh} ATLAS detector-simulation programs~\cite{SOFT-2010-01}, since they were used to unfold the measurements to the particle level as described in Section~\ref{sec:unfolding}. 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 $pp$ interactions per bunch crossing, simulated using the \textsc{Pythia A3} tune~\cite{ATL-PHYS-PUB-2016-017}, and the effect of interactions from bunch crossings before or after the one containing the hard interaction.
 
 
 


\section{Jet reconstruction and event selection}
\label{sec:data_selection}
Jets are reconstructed using the anti-$k_t$ algorithm~\cite{antikt} with radius parameter $R=0.4$ using the \textsc{FastJet} program~\cite{Fastjet}. The inputs to the jet algorithm are particle-flow objects~\cite{pflow}, which make use of both the calorimeter and the inner-detector information to precisely determine the momenta of the input particles. The jet calibration procedure corrects for effects such as energy loss in inactive material, shower leakage, inefficiencies in energy clustering and jet reconstruction. This is performed using a simulation-based correction, in bins of $\eta$ and $\pt$, derived from the relation between the reconstructed jet energy and the energy of the corresponding particle-level jet, not including muons or non-interacting particles. The procedure also includes energy corrections for pile-up, and angular corrections. In a final step, an in situ calibration corrects for residual differences in the jet response and resolution between the MC simulation and the data using $\pt$-balance techniques for dijet, $\gamma+$jet, $\Zboson$+jet and multijet final states~\cite{JES1,JES2}.
 
Events were selected using a single-jet trigger~\cite{TRIG-2016-01} with a minimum $\pt$ threshold of 460~\GeV, to select high-$\pt$ jets while avoiding prescaled triggers. The selected jets are required to have $\pt > 60~\GeV$ and $|\eta| < 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~\cite{ATLAS-CONF-2015-029}. 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, $H_{\mathrm{T}2} = p_{\mathrm{T}1} + p_{\mathrm{T}2} > 1~\TeV$. This ensures a trigger efficiency of $\approx$100\%. About 57.5 million events in data satisfy the selection criteria. The data is binned in ten intervals of $H_{\mathrm{T}2}$ in which the TEEC and ATEEC distributions are measured.


\section{Unfolding to particle level}
\label{sec:unfolding}
To make detector-independent comparisons with particle-level parton shower MC predictions and fixed-order 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 $c\tau>$ 10 mm; these particles are referred to as `stable'. Particle-level jets are reconstructed using the anti-$k_t$ algorithm with $R=0.4$ using stable particles, excluding muons and neutrinos. The fiducial phase space closely follows the event selection criteria defined in Section~\ref{sec:data_selection}. Particle-level jets are required to have $\pt > 60~\GeV$ and $|\eta| < 2.4$. Events with at least two particle-level jets are considered. These events are also required to have $H_{\mathrm{T}2} > 1~\TeV$.
 
The unfolding is performed using an iterative algorithm based on Bayes theorem~\cite{bayes}, 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, $M_{ij}$, obtained from MC simulation, that parameterises the probability of a jet pair generated in bin $i$ being reconstructed in bin $j$. 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. The inverse problem can be written as a linear equation
 
\[
R_i = \sum_{j=1}^N \frac{\mathcal{E}_j}{\mathcal{P}_i} M_{ij} T_j.
\]
 
The quantities $R_i$ and $T_j$ are the contents of the detector-level distribution in bin $i$ and the contents of the particle-level distribution in bin $j$, respectively, while the factors $\mathcal{E}_j$ and $\mathcal{P}_i$ are efficiency and purity corrections, which are estimated by using the MC simulation.
 
The unfolding is performed separately in each bin of $H_{\mathrm{T}2}$ only considering the dependence on $\cos\phi$. The effect of migrations across different $H_{\mathrm{T}2}$ bins is small and is treated by the efficiency and purity corrections. It has been verified that using a two-dimensional transfer matrix, including $H_{\mathrm{T}2}$ and $\cos\phi$, gives results compatible within the statistical uncertainties.
 
Five iterations are necessary for the procedure to converge. 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 \textsc{Pythia 8} are used to obtain the nominal differential cross-sections, whereas \textsc{Sherpa} and \textsc{Herwig 7} MC predictions are used to estimate the systematic uncertainty due to the model dependence, as discussed in Section~\ref{sec:systematics}.
 


\section{Experimental uncertainties}
\label{sec:systematics}
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:
\begin{itemize}
\item{} Jet Energy Scale and Resolution: The jet energy scale (JES) and jet energy resolution (JER) uncertainties are estimated as described in Ref.~\cite{JES1,JES2}. The JES uncertainties consist of a set of 113 independent components, which depend on the jet $\pt$ and $\eta$. 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~\cite{JES1,JES2} are also included. The total JES uncertainty in the $\pt$ value of individual jets varies from approximately 2\% at $\pt$~=~60~\GeV\ to below 1\% for $\pt > 500~\GeV$, with a mild dependence on $\eta$. 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~\cite{JES2}. The effect of the total JER uncertainty is evaluated by smearing the energy of the jets in the MC simulation by about 2\% at $\pt$ = 60~\GeV\ to about 0.5\% for $\pt$ of several hundred \GeV. Each JES and JER uncertainty component is independently propagated by varying the energy and $\pt$ of each jet by one standard deviation. The total uncertainty in the TEEC distribution ranges from a maximum of 2\% for $1~\TeV < H_{\mathrm{T2}} < 1.2~\TeV$ to 1.5\% for $H_{\mathrm{T}2} > 3.5~\TeV$ while, for the ATEEC, it is below 1\% in the whole phase space of the analysis.
 
\item{} Jet Angular Resolution: The jet angular resolution (JAR) uncertainty is conservatively estimated by smearing the azimuthal coordinate $\phi$ of the jets by their resolution in MC simulation. The $\phi$ variations are performed with the $p_{\textrm{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.
 
\item{} 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~\ref{sec:unfolding}. 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.
 
\item{} Monte Carlo Modelling: The modelling uncertainty is estimated by comparing the distributions obtained by using \textsc{Pythia 8}, \textsc{Sherpa} and \textsc{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\phi < -0.95$.
\end{itemize}
 
A breakdown of the relative systematic uncertainties in the measurement of the TEEC and ATEEC distributions is shown in Figure~\ref{fig:syst}.
 
\begin{figure}[H]
\centering
\subfloat[]{\includegraphics[width=0.49\textwidth]{fig_01a.pdf}}
\subfloat[]{\includegraphics[width=0.49\textwidth]{fig_01b.pdf}}
\caption{Breakdown of the systematic uncertainties as a function of (a) the TEEC and (b) the ATEEC for the inclusive $H_{\mathrm{T}2}$ sample, including jet energy scale and resolution, jet angular resolution and uncertainties in the unfolding procedure. Due to the normalisation of the TEEC distribution, the systematic uncertainties show a sharp change in magnitude close to $\cos\phi \simeq 1.0$ for the TEEC and to $\cos\phi\simeq-1.0$ for the ATEEC.}
\label{fig:syst}
\end{figure}
 
 
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~\cite{PERF-2012-01}. 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.


\section{Experimental results}
The unfolded TEEC and ATEEC distributions are presented and compared with the MC predictions described in Section~\ref{sec:MC_simulation}. Figure~\ref{fig:resultsMC1} shows this comparison for the TEEC and ATEEC distributions obtained in the inclusive region $H_{\mathrm{T}2} > 1~\TeV$. Additionally, Figure~\ref{fig:resultsMC2} shows this comparison for both the TEEC and ATEEC for two different bins of $H_{\mathrm{T}2}$.
 
\begin{figure}[H]
\centering
\subfloat[]{\includegraphics[width=0.78\textwidth]{fig_02a.pdf}}\\
 
\subfloat[]{\includegraphics[width=0.78\textwidth]{fig_02b.pdf}}
\caption{Comparison between the data and the MC predictions for (a) the TEEC and (b) the ATEEC for the inclusive $H_{\mathrm{T}2}$ 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.}
\label{fig:resultsMC1}
\end{figure}
 
\begin{figure}[H]
\centering
\subfloat[]{\includegraphics[width=0.44\textwidth]{fig_03a.pdf}}
\subfloat[]{\includegraphics[width=0.44\textwidth]{fig_03b.pdf}}
 
\subfloat[]{\includegraphics[width=0.44\textwidth]{fig_03c.pdf}}
\subfloat[]{\includegraphics[width=0.44\textwidth]{fig_03d.pdf}}
\caption{Comparison between the data and the MC expectations for (a, b) the TEEC and (c, d) the ATEEC for two selected regions of $H_{\mathrm{T}2}$. 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 values of the first bins of the ATEEC distributions are negative and, therefore, not represented in logarithmic scale.}
\label{fig:resultsMC2}
\end{figure}
The TEEC and ATEEC distributions presented in Figures~\ref{fig:resultsMC1} and~\ref{fig:resultsMC2} show the features observed in previous results~\cite{STDM-2014-03, STDM-2016-10}. The TEEC shows two peaks arising from back-to-back ($\cos\phi = -1$) and collinear ($\cos\phi = 1$) configurations, together with a central plateau arising from the wide-angle radiation. The effect of the jet radius $R$ is seen as a kink in the TEEC distributions at $\cos\phi \simeq 0.92$. The ATEEC exhibits a steep fall-off of several orders of magnitude between $\cos\phi = -1$ and $\cos\phi = 0$.
 
At low values of $H_{\mathrm{T}2}$, the best description of the TEEC distributions is achieved by both \textsc{Sherpa} and \textsc{Herwig 7} with the angle-ordered parton shower. The dipole parton shower in \textsc{Herwig 7} does not describe the data well, in particular in the region $|\cos\phi| > 0.4$, where discrepancies of up to 20\% are observed close to the edges of the TEEC distribution. \textsc{Pythia 8} shows discrepancies compared with data for $-0.7 < \cos\phi < -0.2$, where differences of $\mathcal{O}(5\%-10\%)$ are observed. A good description of the ATEEC distribution is achieved by both \textsc{Sherpa} and \textsc{Herwig 7}, while \textsc{Pythia 8} underestimates the measured distribution for $\cos\phi > -0.7$ and the \textsc{Herwig 7} sample matched to the dipole parton shower underestimates the values of the ATEEC for $\cos \phi < -0.7$.
 
For large values of $H_{\mathrm{T}2}$, \textsc{Pythia 8} gives the best description of the TEEC data, while both \textsc{Sherpa} and \textsc{Herwig 7} tend to overestimate the height of the central plateau. The dipole parton shower in \textsc{Herwig 7} shows a similar behaviour to the low $H_{\mathrm{T}2}$ region, with the larger discrepancies observed near the edges of the TEEC distribution, although limited to the region $|\cos\phi| > 0.7$. For the ATEEC, \textsc{Pythia 8} and \textsc{Sherpa} give a good description of the data, while the \textsc{Herwig 7} sample matched to the angle-ordered parton shower overestimates the values of the ATEEC for $\cos\phi < -0.7$. The cancellation of discrepancies in $\phi$-symmetric regions between the data and the dipole shower in \textsc{Herwig 7} brings a reasonable agreement for the ATEEC in the high energy regime.
 
 


\section{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, $\alpha_{\mathrm{s}}(\mu_\mathrm{R})$, as discussed in Ref.~\cite{mitov2022, mitov2023}. The calculations make use of \textsc{OpenLoops2}~\cite{Buccioni:2019sur}, \textsc{FivePointAmplitudes}~\cite{fivepoint} and \textsc{PentagonFunctions++}~\cite{pentagonfunctions} 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 $\Lambda_{\mathrm{QCD}}$ when crossing the corresponding mass thresholds. On the other hand, the perturbative corrections to the three-jet cross-section assume $n_{\mathrm{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.
 
The value of $\alpha_{\mathrm{s}}(\mu_\mathrm{R})$ that enters in the partonic matrix-element calculation is obtained  by evolving the same $\alpha_{\mathrm{s}}(m_\Zboson)$ value as used to determine the PDF set under consideration\@. The NNLO PDF sets \textsc{MMHT2014}~\cite{mmht}, \textsc{NNPDF 3.0}~\cite{nnpdf3} and CT14~\cite{ct14}, as provided in the LHAPDF package~\cite{LHAPDF}, are convolved with the partonic cross-sections.
 
Jets are reconstructed using the anti-$k_{t}$ algorithm with $R=0.4$ as implemented in \textsc{FastJet} \cite{Fastjet}, and the azimuthal angular range in the calculation is restricted to $|\cos\phi| < 0.92$. This helps suppress the sensitivity to resummation effects in the collinear and back-to-back regions.
 
The NNLO corrections are calculated using $\mathcal{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~\cite{NNLOscheme1,NNLOscheme2}. For illustrative purposes, Figure~\ref{fig:nnloDiagrams} shows sample diagrams corresponding to real--real, real--virtual and virtual--virtual contributions to the $gg\to \mathrm{jjj}$ process at NNLO.
 
\begin{figure}[H]
\centering
\subfloat[]{\includegraphics[width=0.23\textwidth]{fig_04a.pdf}}
\hspace{1.2cm}
\subfloat[]{\includegraphics[width=0.23\textwidth]{fig_04b.pdf}}
\hspace{1.2cm}
\subfloat[]{\includegraphics[width=0.23\textwidth]{fig_04c.pdf}}
\caption{Feynman diagrams contributing to the $gg\to \mathrm{jjj}$ process at NNLO accuracy in pQCD\@. These particular configurations contribute to the (a) real--real, (b) real--virtual, and (c) virtual--virtual finite terms.}
\label{fig:nnloDiagrams}
\end{figure}
For each event, the renormalisation and factorisation scales are set to the scalar sum of the transverse momenta of all final-state partons, $\mu_R = \mu_F = \hat{H}_{\mathrm{T}}$~\cite{CERN-TH-2018-159, mitov2022, mitov2023}. This choice allows the effect of the third jet to be taken into account in the definition of the scale, which was not included in previous analyses~\cite{STDM-2014-03, STDM-2016-10}, where $\mu_R = \mu_F = H_{\mathrm{T}2}/2$ was used.
 
The pQCD predictions for the TEEC and ATEEC functions are calculated at parton level only. 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 \textsc{Pythia~8}~\cite{ATL-PHYS-PUB-2014-021}, and are close to unity to within 0.5\%, except for $\cos\phi > 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.
\begin{itemize}
\item{} The scale uncertainties are evaluated by considering independent variations of the renormalisation and factorisation scales by a factor of two, excluding those in which $\mu_\mathrm{R}$ and $\mu_\mathrm{F}$ are varied in opposite directions. 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. This is understood to be because at NNLO, two different terms of the form $c_k\log^k(\mu_R/Q)$, with $k=1,2$, contribute to the cross-section dependence on $\mu_R$. Thus, the asymmetry in the uncertainty band depends on the interplay between the coefficients $c_k$ of both terms. 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.
\item{} The PDF uncertainties are obtained by considering the set of eigenvectors / replicas provided by each PDF group~\cite{mmht,nnpdf3,ct14}, 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 $H_{\mathrm{T2}}$ bins, and can reach up to 2\% for the higher $H_{\mathrm{T2}}$ bins.
\item{} The uncertainties in the non-perturbative corrections are evaluated by considering the envelope of the differences between five different predictions and the \textsc{Pythia 8} A14 tune~\cite{ATL-PHYS-PUB-2014-021}. These predictions were generated using \textsc{Pythia 8.235} and \textsc{Herwig 7.2.1} with the angle-ordered parton shower. They include the \textsc{4C}~\cite{Tune4C} and \textsc{Monash}~\cite{Skands:2014pea} \textsc{Pythia 8} tunes, and the \textsc{Herwig 7} default, soft and baryonic reconnection tunes~\cite{herwigTunes}. The value of this uncertainty is below 1\% in the phase-space region $|\cos\phi| < 0.92$, where $\alpha_{\mathrm{s}}$ is determined.
\end{itemize}
The relative values of the theoretical uncertainties for the TEEC and ATEEC distributions are shown in Figure~\ref{fig:theoSys} for the inclusive $H_{\mathrm{T2}}$ bin. The effect of varying $\alpha_s(m_\Zboson)$ within the range 0.117 -- 0.119 is also shown for illustrative purposes.
\begin{figure}[H]
\subfloat[]{\includegraphics[width=0.49\textwidth]{fig_05a.pdf}}
\subfloat[]{\includegraphics[width=0.49\textwidth]{fig_05b.pdf}}
\caption{Theoretical uncertainties in the NNLO pQCD prediction for (a) the TEEC and (b) the ATEEC distributions for $H_{\mathrm{T}2} > 1$~\TeV. The predictions use the \textsc{MMHT2014} PDF\@. The effect of a variation of $\alpha_s(m_\Zboson)$ from 0.118 by $\pm 0.001$ is also shown for illustration.}
\label{fig:theoSys}
\end{figure}
Figures~\ref{fig:ratioTEEC} and~\ref{fig:ratioATEEC} 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\phi > 0.80$, shows a large improvement for all $H_{\mathrm{T2}}$ bins. For the higher $H_{\mathrm{T2}}$ bins, and for the \textsc{MMHT2014} PDF, the NNLO theoretical predictions are slightly above the data. Since the \textsc{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 $x$. The results obtained using the \textsc{CT18} PDF set agree well with those for \textsc{MMHT2014}.
 
\begin{figure}[H]
\vspace{1.8cm}
\centering
\includegraphics[width=0.955\textwidth]{fig_06.pdf}
\caption{Ratios of the theoretical predictions for the TEEC functions at LO and NNLO to the NLO calculations, together with the ratios of the data to NLO predictions. The hatched band, where visible, shows the statistical uncertainty in the NNLO prediction. The predictions use the \textsc{MMHT2014} PDF, where the value of the strong coupling constant is set to $\alpha_{\mathrm{s}}(m_\Zboson) = 0.1180$. The uncertainty bands correspond to the scale uncertainties for each perturbative order.}
\label{fig:ratioTEEC}
\end{figure}
\begin{figure}[H]
\vspace{1.8cm}
\centering
\includegraphics[width=0.955\textwidth]{fig_07.pdf}
\caption{Ratios of the theoretical predictions for the ATEEC functions at LO and NNLO to the NLO calculations, together with the ratios of the data to NLO predictions. The hatched band, where visible, shows the statistical uncertainty in the NNLO prediction. The predictions use the \textsc{MMHT2014} PDF, where the value of the strong coupling constant is set to $\alpha_{\mathrm{s}}(m_\Zboson) = 0.1180$. The uncertainty bands correspond to the scale uncertainties for each perturbative order.}
\label{fig:ratioATEEC}
\end{figure}
\clearpage
 


\section{Determination of the strong coupling constant}
The strong coupling constant at the scale of the pole mass of the $\Zboson$ boson, $\alpha_{\mathrm{s}}(m_\Zboson)$ can be determined from the comparison of the data with the theoretical predictions by considering the following $\chi^2$ function
\begin{equation}
\chi^2(\alpha_{\mathrm{s}},\vec{\lambda}) = \sum_{\mathrm{bins}}\frac{(x_i-F_i(\alpha_{\mathrm{s}},\vec{\lambda}))^2}{\Delta x_i^2 + \Delta \xi_i^2} + \sum_{k}\lambda_k^2,
\label{eq:chi2a}
\end{equation}
where the theoretical predictions are varied according to
\begin{equation}
F_{i}(\alpha_{\mathrm{s}},\vec{\lambda}) = \psi_{i}(\alpha_{\mathrm{s}})\left(1+\sum_k\lambda_k\sigma_k^{(i)}\right).
\label{eq:chi2b}
\end{equation}
In Eq.~(\ref{eq:chi2a}) and Eq.~(\ref{eq:chi2b}), $\alpha_{\mathrm{s}}$ stands for $\alpha_{\mathrm{s}}(m_\Zboson)$; $x_i$ is the value of the $i$-th point of the distribution as measured in data, while $\Delta x_i$ is its statistical uncertainty. The statistical uncertainty in the theoretical predictions dominates over the statistical uncertainty in data, and is also included as $\Delta \xi_i$, while $\sigma_k^{(i)}$ is the relative value of the $k$-th correlated source of systematic uncertainty in bin $i$.
 
This technique takes into account the correlations between the different sources of systematic uncertainty discussed in Section~\ref{sec:systematics} by introducing the nuisance parameters $\{\lambda_k\}$, one for each source of experimental uncertainty, each of which follow a standard normal distribution. Thus, the minimum of the $\chi^2$ function defined in Eq.~(\ref{eq:chi2a}) is found in a 150-dimensional space, in which 149 correspond to nuisance parameters $\left\{\lambda_i\right\}$ and one to $\alpha_{\mathrm{s}}(m_\Zboson)$. 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 \textsc{Herwig 7} and \textsc{Sherpa}. The envelope of the $\alpha_s$ 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 $\psi_{i}(\alpha_{\mathrm{s}})$ for bin $i$. For each PDF set, the corresponding $\alpha_{\mathrm{s}}(m_\Zboson)$ variation range is considered and the theoretical prediction is obtained for each value of $\alpha_{\mathrm{s}}(m_\Zboson)$. The functions $\psi_{i}(\alpha_{\mathrm{s}})$ are then obtained by fitting the predicted values of the TEEC (ATEEC) in each $(H_{\mathrm{T}2}, \cos\phi)$ bin to a third-order polynomial in $\alpha_s$.
 
For both the TEEC and ATEEC functions, the fits to extract $\alpha_{\mathrm{s}}(m_\Zboson)$ are repeated separately for each $H_{\mathrm{T}2}$ interval, thus determining a value of $\alpha_{\mathrm{s}}(m_\Zboson)$ for each energy bin. The theoretical uncertainties are determined by shifting the theory distributions by each of the uncertainties separately, recalculating the functions $\psi_i(\alpha_{\mathrm{s}})$ and determining a new value of $\alpha_{\mathrm{s}}(m_\Zboson)$. 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 $\alpha_{\mathrm{s}}(m_\Zboson)$, which are obtained in the $\overline{\mbox{MS}}$ subtraction scheme, is then evolved to the corresponding measured scale using the NNLO solution to the RGE\@. When evolving $\alpha_{\mathrm{s}}(m_\Zboson)$ to $\alpha_{\mathrm{s}}(Q)$, the appropriate matching conditions for the strong coupling constant at the $n_{\mathrm{f}}=5$ and $n_{\mathrm{f}}=6$ flavour thresholds are applied so that $\alpha_{\mathrm{s}}(Q)$ is a continuous function across quark thresholds. The value of $Q$ is obtained as the average value of $\mu_R = \hat{H}_\mathrm{T}$ calculated at NNLO for each $H_{\mathrm{T}2}$ bin.
 
 
The values of $\alpha_s(m_\Zboson)$ obtained from global fits to the TEEC and ATEEC distributions, using the \textsc{MMHT 2014}, \textsc{CT14} and \textsc{NNPDF 3.0} PDF sets are shown in Tables~\ref{tab:aspdfTEEC} and~\ref{tab:aspdfATEEC}, together with the corresponding $\chi^2$ values. The values of $\alpha_{\mathrm{s}}(m_\Zboson)$ determined from fits to the TEEC functions, both in the inclusive and exclusive $H_{\mathrm{T}2}$ bins, are presented in Table~\ref{tab:teecAlpha}, while the values obtained from fits to the ATEEC functions are presented in Table~\ref{tab:ateecAlpha}. The average $Q$ values and the reduced $\chi^2$ values at the minimum are also shown. The values displayed in Tables~\ref{tab:teecAlpha} and~\ref{tab:ateecAlpha} are obtained using the \textsc{MMHT 2014} PDF in the theoretical predictions.
 
\begin{table}[H]
\caption{Values of $\alpha_s(m_\Zboson)$ obtained from fits to the TEEC distributions to the NNLO theoretical predictions obtained using different PDF sets. The uncertainty referred to as $\mu$ 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 $\chi^2$ function at its minimum, as well as the number of degrees of freedom is also indicated.}
\label{tab:aspdfTEEC}
\centering
\footnotesize
\vspace{0.5cm}
\begin{tabular}{ c c c }
\hline
\T\B PDF & $\alpha_s(m_\Zboson)$ value & $\chi^2 / N_\mathrm{dof}$\\
\hline
\T MMHT 2014 & $0.1175 \pm 0.0001 \mbox{ (stat.)} \pm 0.0006 \mbox{ (sys.)}^{+0.0032}_{-0.0011} \mbox{ ($\mu$)} \pm 0.0011 \mbox{ (PDF)} \pm 0.0002 \mbox{ (NP)} \pm 0.0005 \mbox{ (mod.)}$ & 318 / 251\\
\T CT14 & $0.1196 \pm 0.0001 \mbox{ (stat.)} \pm 0.0006 \mbox{ (sys.)}^{+0.0035}_{-0.0010} \mbox{ ($\mu$)} \pm 0.0016 \mbox{ (PDF)} \pm 0.0002 \mbox{ (NP)} \pm 0.0006 \mbox{ (mod.)}$ & 262 / 251\\
\T\B NNPDF 3.0 & $0.1191 \pm 0.0001 \mbox{ (stat.)} \pm 0.0006 \mbox{ (sys.)}^{+0.0040}_{-0.0011} \mbox{ ($\mu$)} \pm 0.0020 \mbox{ (PDF)} \pm 0.0003 \mbox{ (NP)} \pm 0.0007 \mbox{(mod.)}$ & 300 / 251\\
\hline
\end{tabular}
\end{table}
 
\begin{table}[H]
\caption{Values of $\alpha_s(m_\Zboson)$ obtained from fits to the ATEEC distributions to the NNLO theoretical predictions obtained using different PDF sets. The uncertainty referred to as $\mu$ 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 $\chi^2$ function at its minimum, as well as the number of degrees of freedom is also indicated.}
\label{tab:aspdfATEEC}
\centering
\footnotesize
\vspace{0.5cm}
\begin{tabular}{ c c c }
\hline
\T\B PDF & $\alpha_s(m_\Zboson)$ value & $\chi^2 / N_\mathrm{dof}$\\
\hline
\T MMHT 2014 & $0.1185 \pm 0.0005 \mbox{ (stat.)} \pm 0.0008 \mbox{ (sys.)}^{+0.0022}_{-0.0002} \mbox{ ($\mu$)} \pm 0.0011 \mbox{ (PDF)} \pm 0.0004 \mbox{ (NP)} \pm 0.0001 \mbox{ (mod.)}$ & 110 / 117\\
\T CT14 & $0.1200 \pm 0.0006 \mbox{ (stat.)} \pm 0.0009 \mbox{ (sys.)}^{+0.0027}_{-0.0001} \mbox{ ($\mu$)} \pm 0.0016 \mbox{ (PDF)} \pm 0.0005 \mbox{ (NP)} \pm 0.0001 \mbox{ (mod.)}$ & 110 / 117\\
\T\B NNPDF 3.0 & $0.1199 \pm 0.0006 \pm \mbox{ (stat.)} 0.0009 \mbox{ (sys.)}^{+0.0027}_{-0.0002} \mbox{ ($\mu$)} \pm 0.0017 \mbox{ (PDF)} \pm 0.0005 \mbox{ (NP)} \pm 0.0001 \mbox{(mod.)}$ & 108 / 117\\
\hline
\end{tabular}
\end{table}
 
The final values of $\alpha_s$ obtained from global fits to the TEEC and ATEEC distributions using the \textsc{MMHT 2014} PDF are
\[
\begin{array}{c}
\B\alpha_\mathrm{s}(m_\Zboson) = 0.1175 \pm 0.0006 \mbox{ (exp.)} ^{+0.0034}_{-0.0017} \mbox{ (theo.) and}\\
\T\alpha_\mathrm{s}(m_\Zboson) = 0.1185 \pm 0.0009 \mbox{ (exp.)} ^{+0.0025}_{-0.0012} \mbox{ (theo.)}.
\end{array}
\]
The central values of $\alpha_s$ determined from both observables are correlated with a Pearson correlation coefficient $\rho=0.86\pm 0.02 \mbox{ (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 $\alpha_{\mathrm{s}}(m_\Zboson)$ 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~\ref{fig:teecRatio} 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~\ref{fig:ateecRatio} 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 $\alpha_{\mathrm{s}}(m_\Zboson)$ obtained in Tables~\ref{tab:teecAlpha} and~\ref{tab:ateecAlpha} 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~\ref{fig:teecRunning} and~\ref{fig:ateecRunning} show the values of $\alpha_{\mathrm{s}}(Q)$ together with the world average band provided by the Particle Data Group~\cite{PDG2022} and values of $\alpha_\mathrm{s}$ obtained in other analyses~\cite{D0-1, D0-2, STDM-2014-03, STDM-2016-10, D0-1, STDM-2012-19, CMS-SMP-14-001, CMS-SMP-12-027, CMS-QCD-11-003, CMS-TOP-18-004}. To compare the results with other experiments, the value of $Q$ is chosen as half of the average value of $\hat{H}_\mathrm{T}$ for each $H_{\mathrm{T}2}$ bin. The results show a good agreement between all measurements and the RGE prediction. The values of $\alpha_s$ obtained from fits to the theoretical predictions using different PDF sets, shown in Figs.~\ref{fig:teecRunningPDF} and~\ref{fig:ateecRunningPDF}, are found to be in good agreement. Since the main purpose of the present analysis is to provide a first determination of $\alpha_s$ using NNLO predictions for three-jet cross sections, a combined determination of $\alpha_s$ and the PDF~\cite{forte} is not attempted. Such analysis would require to use very large computing resources and is well out of the scope of this paper.
 
\newpage
\begin{landscape}
\begin{table}[H]
\caption{Values of the strong coupling constant at the $Z$ boson mass scale, $\alpha_{\mathrm{s}}(m_\Zboson)$, obtained from fits to the TEEC function using \textsc{MMHT 2014} parton distribution functions. The value referred to as `inclusive' is obtained from a fit to the TEEC distribution for $H_{\mathrm{T}2} > 1~\TeV$. The values of the average interaction scale $\langle Q \rangle$ are shown in the first column, while the values of the $\chi^{2}$ function at the minimum are shown in the third column. The uncertainty referred to as $\mu$ 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.}
\label{tab:teecAlpha}
\small
\vspace{2.0cm}
\centering
\begin{tabular}{ c c c }
\hline
\T\B $\langle \hat{H}_{\mathrm{T}} \rangle$ [\GeV] & $\alpha_{\mathrm{s}}(m_\Zboson)$ value (\textsc{MMHT 2014}) & $\chi^{2} / N_{\mathrm{dof}}$ \\
\hline
\T\B Inclusive & 0.1188 $\pm$ 0.0002 (stat.) $\pm$ 0.0007 (syst.) $^{+0.0030}_{-0.0002}$ ($\mu$) $\pm$ 0.0010 (PDF) $\pm$ 0.0002 (NP) $\pm$ 0.0008 (mod.) & 16.0 / 27 \\
\hline
\T\B 1302 & 0.1186 $\pm$ 0.0003 (stat.) $\pm$ 0.0009 (syst.) $^{+0.0031}_{-0.0001}$ ($\mu$) $\pm$ 0.0010 (PDF) $\pm$ 0.0003 (NP) $\pm$ 0.0007 (mod.) & 16.4 / 27 \\
\T\B 1518 & 0.1182 $\pm$ 0.0003 (stat.) $\pm$ 0.0009 (syst.) $^{+0.0027}_{-0.0005}$ ($\mu$) $\pm$ 0.0010 (PDF) $\pm$ 0.0003 (NP) $\pm$ 0.0008 (mod.) & 15.1 / 27 \\
\T\B 1732 & 0.1191 $\pm$ 0.0003 (stat.) $\pm$ 0.0011 (syst.) $^{+0.0030}_{-0.0005}$ ($\mu$) $\pm$ 0.0011 (PDF) $\pm$ 0.0005 (NP) $\pm$ 0.0010 (mod.) & 19.8 / 27 \\
\T\B 1944 & 0.1179 $\pm$ 0.0003 (stat.) $\pm$ 0.0011 (syst.) $^{+0.0030}_{-0.0005}$ ($\mu$) $\pm$ 0.0011 (PDF) $\pm$ 0.0005 (NP) $\pm$ 0.0009 (mod.) & 21.3 / 27 \\
\T\B 2153 & 0.1175 $\pm$ 0.0004 (stat.) $\pm$ 0.0012 (syst.) $^{+0.0029}_{-0.0004}$ ($\mu$) $\pm$ 0.0011 (PDF) $\pm$ 0.0003 (NP) $\pm$ 0.0009 (mod.) & 32.9 / 27 \\
\T\B 2396 & 0.1179 $\pm$ 0.0003 (stat.) $\pm$ 0.0012 (syst.) $^{+0.0029}_{-0.0005}$ ($\mu$) $\pm$ 0.0012 (PDF) $\pm$ 0.0011 (NP) $\pm$ 0.0012 (mod.) & 27.5 / 27 \\
\T\B 2706 & 0.1164 $\pm$ 0.0004 (stat.) $\pm$ 0.0015 (syst.) $^{+0.0030}_{-0.0005}$ ($\mu$) $\pm$ 0.0013 (PDF) $\pm$ 0.0005 (NP) $\pm$ 0.0011 (mod.) & 35.1 / 27 \\
\T\B 3042 & 0.1162 $\pm$ 0.0005 (stat.) $\pm$ 0.0017 (syst.) $^{+0.0031}_{-0.0005}$ ($\mu$) $\pm$ 0.0013 (PDF) $\pm$ 0.0002 (NP) $\pm$ 0.0015 (mod.) & 33.1 / 27 \\
\T\B 3476 & 0.1141 $\pm$ 0.0007 (stat.) $\pm$ 0.0017 (syst.) $^{+0.0033}_{-0.0011}$ ($\mu$) $\pm$ 0.0014 (PDF) $\pm$ 0.0002 (NP) $\pm$ 0.0020 (mod.) & 15.1 / 13 \\
\T\B 4189 & 0.1116 $\pm$ 0.0011 (stat.) $\pm$ 0.0018 (syst.) $^{+0.0030}_{-0.0009}$ ($\mu$) $\pm$ 0.0015 (PDF) $\pm$ 0.0002 (NP) $\pm$ 0.0020 (mod.) & 14.0 / 13 \\
\hline
\end{tabular}
\end{table}
\end{landscape}
 
\begin{landscape}
\begin{table}[H]
\caption{Values of the strong coupling constant at the $Z$ boson mass scale, $\alpha_{\mathrm{s}}(m_\Zboson)$, obtained from fits to the ATEEC function using \textsc{MMHT 2014} parton distribution functions. The value referred to as `inclusive' is obtained from a fit to the ATEEC distribution for $H_{\mathrm{T}2} > 1~\TeV$. The values of the average interaction scale $\langle Q \rangle$ are shown in the first column, while the values of the $\chi^{2}$ function at the minimum are shown in the third column. The uncertainty referred to as $\mu$ 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.}
\label{tab:ateecAlpha}
\small
\vspace{2.0cm}
\centering
\begin{tabular}{ c c c }
\hline
\T\B $\langle \hat{H}_{\mathrm{T}} \rangle$ [\GeV] & $\alpha_{\mathrm{s}}(m_\Zboson)$ value (\textsc{MMHT 2014}) & $\chi^{2} / N_{\mathrm{dof}}$ \\
\hline
\T\B Inclusive & 0.1194 $\pm$ 0.0009 (stat.) $\pm$ 0.0007 (syst.) $^{+0.0023}_{-0.0000}$ ($\mu$) $\pm$ 0.0011 (PDF) $\pm$ 0.0005 (NP) $\pm$ 0.0000 (mod.) & 10.7 / 12 \\
\hline
\T\B 1302 & 0.1195 $\pm$ 0.0011 (stat.) $\pm$ 0.0006 (syst.) $^{+0.0024}_{-0.0000}$ ($\mu$) $\pm$ 0.0011 (PDF) $\pm$ 0.0006 (NP) $\pm$ 0.0000 (mod.) & 10.6 / 12 \\
\T\B 1518 & 0.1191 $\pm$ 0.0011 (stat.) $\pm$ 0.0007 (syst.) $^{+0.0020}_{-0.0001}$ ($\mu$) $\pm$ 0.0011 (PDF) $\pm$ 0.0004 (NP) $\pm$ 0.0001 (mod.) & 7.2 / 12 \\
\T\B 1732 & 0.1187 $\pm$ 0.0015 (stat.) $\pm$ 0.0009 (syst.) $^{+0.0026}_{-0.0003}$ ($\mu$) $\pm$ 0.0012 (PDF) $\pm$ 0.0010 (NP) $\pm$ 0.0002 (mod.) & 7.2 / 12 \\
\T\B 1944 & 0.1178 $\pm$ 0.0016 (stat.) $\pm$ 0.0009 (syst.) $^{+0.0022}_{-0.0001}$ ($\mu$) $\pm$ 0.0013 (PDF) $\pm$ 0.0007 (NP) $\pm$ 0.0000 (mod.) & 11.2 / 12 \\
\T\B 2153 & 0.1174 $\pm$ 0.0017 (stat.) $\pm$ 0.0009 (syst.) $^{+0.0022}_{-0.0002}$ ($\mu$) $\pm$ 0.0013 (PDF) $\pm$ 0.0007 (NP) $\pm$ 0.0001 (mod.) & 12.8 / 12 \\
\T\B 2396 & 0.1187 $\pm$ 0.0017 (stat.) $\pm$ 0.0010 (syst.) $^{+0.0017}_{-0.0000}$ ($\mu$) $\pm$ 0.0012 (PDF) $\pm$ 0.0007 (NP) $\pm$ 0.0004 (mod.) & 11.7 / 12 \\
\T\B 2706 & 0.1148 $\pm$ 0.0026 (stat.) $\pm$ 0.0014 (syst.) $^{+0.0024}_{-0.0001}$ ($\mu$) $\pm$ 0.0015 (PDF) $\pm$ 0.0007 (NP) $\pm$ 0.0002 (mod.) & 18.3 / 12 \\
\T\B 3042 & 0.1169 $\pm$ 0.0031 (stat.) $\pm$ 0.0012 (syst.) $^{+0.0018}_{-0.0000}$ ($\mu$) $\pm$ 0.0015 (PDF) $\pm$ 0.0014 (NP) $\pm$ 0.0009 (mod.) & 13.0 / 12 \\
\T\B 3476 & 0.1141 $\pm$ 0.0052 (stat.) $\pm$ 0.0016 (syst.) $^{+0.0025}_{-0.0007}$ ($\mu$) $\pm$ 0.0018 (PDF) $\pm$ 0.0008 (NP) $\pm$ 0.0011 (mod.) & 4.9 / 6 \\
\T\B 4189 & 0.1096 $\pm$ 0.0085 (stat.) $\pm$ 0.0009 (syst.) $^{+0.0013}_{-0.0000}$ ($\mu$) $\pm$ 0.0009 (PDF) $\pm$ 0.0002 (NP) $\pm$ 0.0007 (mod.) & 6.2 / 6 \\
\hline
\end{tabular}
\end{table}
\end{landscape}
 
\clearpage
\begin{figure}[H]
\vspace{2.cm}
\centering
\includegraphics[width=1.00\textwidth]{fig_08.pdf}
\caption{Ratios of the data to the fitted theoretical predictions for the TEEC measurements in inclusive and exclusive $H_{\mathrm{T}2}$ bins. The solid band shows the theoretical uncertainties, dominated by the scale variations, while the error bars show the experimental uncertainties, where correlations between the fit parameters are taken into account. The hatched band shows the theoretical uncertainties, other than the scale variations. The central values of $\alpha_s(m_Z)$ determined for each $H_{\mathrm{T}2}$ bin are indicated.}
\label{fig:teecRatio}
\end{figure}
 
 
\clearpage
\begin{figure}[H]
\vspace{2.cm}
\centering
\includegraphics[width=1.00\textwidth]{fig_09.pdf}
\caption{Ratios of the data to the fitted theoretical predictions for the ATEEC measurements in inclusive and exclusive $H_{\mathrm{T}2}$ bins. The solid band shows the theoretical uncertainties, dominated by the scale variations, while the error bars show the experimental uncertainties, where correlations between the fit parameters are taken into account. The hatched band shows the theoretical uncertainties, other than the scale variations. The central values of $\alpha_s(m_Z)$ determined for each $H_{\mathrm{T}2}$ bin are indicated.}
\label{fig:ateecRatio}
\end{figure}
 
\clearpage
\begin{figure}[H]
\vspace{-0.3cm}
\centering
\includegraphics[width=0.83\textwidth]{fig_10.pdf}
\caption{Comparison of the values of $\alpha_\mathrm{s}(Q)$ determined from fits to the TEEC functions with the QCD prediction using the world average as input (hatched band) and the value obtained from the global fit (solid band). Results from previous analyses, both from ATLAS and from other experiments, are also included, showing an excellent agreement with the current measurements and with the world average. The value of $Q$ for the TEEC 13~\TeV\ points is chosen as half of the average of $\hat{H}_{\mathrm{T}}$ in each $H_{\mathrm{T}2}$ bin.}
\label{fig:teecRunning}
\end{figure}
 
\begin{figure}[H]
\vspace{-0.3cm}
\centering
\includegraphics[width=0.83\textwidth]{fig_11.pdf}
\caption{Comparison of the values of $\alpha_\mathrm{s}(Q)$ determined from fits to the ATEEC functions with the QCD prediction using the world average as input (hatched band) and the value obtained from the global fit (solid band). Results from previous analyses, both from ATLAS and from other experiments, are also included, showing an excellent agreement with the current measurements and with the world average. The value of $Q$ for the TEEC 13~\TeV\ points is chosen as half of the average of $\hat{H}_{\mathrm{T}}$ in each $H_{\mathrm{T}2}$ bin.}
\label{fig:ateecRunning}
\end{figure}
 
\clearpage
 
\begin{figure}[H]
\vspace{0.7cm}
\centering
\includegraphics[width=0.83\textwidth]{fig_12.pdf}
\caption{Values of $\alpha_s(\mu_R)$ determined in each bin of $H_{\mathrm{T}2}$ from fits to the TEEC distributions using theoretical predictions with the \textsc{MMHT2014}, \textsc{CT14} and \textsc{NNPDF 3.0} PDF sets.}
\label{fig:teecRunningPDF}
\end{figure}
 
\begin{figure}[H]
\vspace{-0.3cm}
\centering
\includegraphics[width=0.83\textwidth]{fig_13.pdf}
\caption{Values of $\alpha_s(\mu_R)$ determined in each bin of $H_{\mathrm{T}2}$ from fits to the ATEEC distributions using theoretical predictions with the \textsc{MMHT2014}, \textsc{CT14} and \textsc{NNPDF 3.0} PDF sets.}
\label{fig:ateecRunningPDF}
\end{figure}
 

\FloatBarrier

\section{Conclusions}
A measurement of transverse energy--energy correlations and their corresponding asymmetries in multijet events using 139~\ifb of $pp$ collisions at $\sqrt{s} = 13~\TeV$ collected with the ATLAS detector at the LHC is presented. High-energy multijet events are selected by requiring the scalar sum of $\pt$ of the two leading jets, $H_{\mathrm{T}2}$, 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 \textsc{Pythia~8}, \textsc{Sherpa} and \textsc{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. From this comparison, values of $\alpha_{\mathrm{s}}(m_\Zboson)$ are obtained for both the $H_{\mathrm{T}2}$-inclusive measurement and for the different exclusive bins in $H_{\mathrm{T}2}$ using a $\chi^2$ fit. The values of $\alpha_\mathrm{s}(m_\Zboson)$ obtained from the global fit to all $H_{\mathrm{T}2}$ bins simultaneously, are
\[
\begin{array}{c}
\B\alpha_\mathrm{s}(m_\Zboson) = 0.1175 \pm 0.0006 \mbox{ (exp.)} ^{+0.0034}_{-0.0017} \mbox{ (theo.) and}\\
\T\alpha_\mathrm{s}(m_\Zboson) = 0.1185 \pm 0.0009 \mbox{ (exp.)} ^{+0.0025}_{-0.0012} \mbox{ (theo.)}
\end{array}
\]
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 $\alpha_s$.
 
The evolution of the values obtained from each of the exclusive fits to the corresponding measurement scales leads to values of $\alpha_\mathrm{s}(Q)$ that are compared with the next-to-next-to-leading-order solution of the RGE. The results show a good agreement with the RGE predictions up to the highest energy scales, and with previous measurements by ATLAS and by other experiments.



 
 
 
\clearpage
\printbibliography
 

 
\clearpage
 

 
\end{document}
