Transverse Energy-Energy Correlations of jets in the electron-proton Deep Inelastic Scattering at HERA

We study the event shape variables, transverse energy energy correlation TEEC $(\cos \phi)$ and its asymmetry ATEEC $(\cos \phi)$ in deep inelastic scattering (DIS) at the electron-proton collider HERA, where $\phi$ is the angle between two jets defined using a transverse-momentum $(k_T)$ jet algorithm. At HERA, jets are defined in the Breit frame, and the leading nontrivial transverse energy energy correlations arise from the 3-jet configurations. With the help of the NLOJET++, these functions are calculated in the leading order (LO) and the next-to-leading order (NLO) approximations in QCD at the electron-proton center-of-mass energy $\sqrt{s}=314$ GeV. We restrict the angular region to $-0.8 \leq \cos \phi \leq 0.8$, as the forward- and backward-angular regions require resummed logarithmic corrections, which we have neglected in this work. Following experimental jet-analysis at HERA, we restrict the DIS-variables $x$, $y=Q^2/(x s)$, where $Q^2=-q^2$ is the negative of the momentum transfer squared $q^2$, to $0 \leq x \leq 1$, $0.2 \leq y \leq 0.6$, and the pseudo-rapidity variable in the laboratory frame $(\eta^{\rm {lab}})$ to the range $-1 \leq \eta^{\rm {lab}} \leq 2.5$. The TEEC and ATEEC functions are worked out for two ranges in $Q^2$, defined by $5.5~{\rm GeV}^2 \leq Q^2 \leq 80~{\rm GeV}^2$, called the low-$Q^2$-range, and $150~{\rm GeV}^2 \leq Q^2 \leq 1000~{\rm GeV}^2$, called the high-$Q^2$-range. We show the sensitivity of these functions on the parton distribution functions (PDFs), the factorization $(\mu_F)$ and renormalization $(\mu_R)$ scales, and on $\alpha_s(M_Z)$. Of these the correlations are stable against varying the scale $\mu_F$ and the PDFs, but they do depend on $\mu_R$. These studies are useful in the analysis of the HERA data, including the determination of $\alpha_s(M_Z)$ from the shape variables.


I. INTRODUCTION
Event shape variables involving the energy-momentum variables of hadrons and jets have played a crucial role in testing Quantum Chromodyamics (QCD), providing a detailed comparison with the experimentally measured shapes in high energy collisions and in determining the strong interaction coupling constant α s (Q 2 ). Of these, the energyenergy correlation (EEC) and its asymmetry (AEEC), introduced by Basham et al. in e + e − annihilation [1,2] have received a lot of experimental and theoretical attention. Next-to-leading order (NLO) corrections in α s (Q 2 ) were calculated long ago for the EEC in e + e − annihilation, using a number of different methods to regulate the soft and collinear divergences [3][4][5][6][7][8][9][10]. Accurate numerical results for the EEC are available from the program Event 2, based on the dipole subtraction technique [11,12]. EEC has also been calculated to NNLO accuracy in perturbative QCD [13,14]. Recent advances in theoretical calculational techniques have led to a renaissance of interest in this topic. In particular, an analytic NLO calculations of the EEC in e + e − annihilation [15,16], and an all-order factorization formula for the EEC in the back-to-back limit [17][18][19][20], are now available. We also mention here the derivation of the EEC function in the maximally supersymmetric N = 4 super-Yang-Mills theory in the NLO accuracy [21], which has been recently extended up to NNLO accuracy [22]. Experimental measurements of EEC in e + e − annihilation are discussed in [23][24][25][26][27].
Following EEC in e + e − annihilation, transverse energy-energy correlation (TEEC) and the corresponding asymmetry (ATEEC) were introduced in hadronic collisions at the SPPS [28], but did not evoke much experimental interest. With the advent of the LHC era, NLO corrections were calculated in pp collisions [29]. They have been used by the ATLAS collaboration for comparison with data and in the determination of α s (M Z ) from these shape functions [30,31]. Recently, TEEC in the dijet back-to-back limit in hadronic collisions has been derived, achieving an impressive perturbative simplicity [32]. Currently the TEEC-data in pp collisions are restricted in their theoretical interpretation to NLO accuracy.
What concerns deep inelastic scattering (DIS), event shape variables have also received a lot of theoretical attention [33][34][35][36][37][38]. Prominent among them are the thrust-distribution, 1-jettiness, jet-broadening, and the C parameter, which have been calculated to very high accuracy in fixed order (NNLO) [38], and in the resummed leading logarithms (N 3 LL) [37]. Some of these event shape variables have been measured by the H1 [39] and ZEUS [40] collaborations at HERA. The definitions of these shape variables together with some others, such as the jet shape, can be seen in [41], where DIS and photoproduction experiments at HERA are reviewed. However, to the best of our knowledge, the transverse energy-energy correlation between the final state jets in deep elastic scattering has neither been calculated nor measured so far. Analogous to the TEEC for hadronic collisions [28,29], TEEC in DIS is introduced in Eq. (1) in the next section. It involves transverse energy correlations in two jets, defined by a jet-definition and jet algorithm, separated by an azimuthal angle φ. We calculate TEEC and its asymmetry in DIS at HERA under realistic experimental conditions. Jets at HERA are defined in the Breit frame, in which the exchanged photon is at rest and the incoming and outgoing quarks are along the z direction. In this frame, the involved hadronic final states have zero total transverse momentum, and thus the leading nontrival transverse energy-energy correlation comes from the 3-jet configurations. To match the measurements of jets at HERA, we adopt the transverse-momentum (k T ) algorithm to classify the jets [42] and calculate the TEEC and its asymmetry (ATEEC) in the kinematic conditions employed typically in H1 and ZEUS. The calculations are done in the NLO accuracy in the central angular region, −0.8 ≤ cos φ ≤ 0.8. This avoids the back-to-back angular configuration, i.e., near φ = π, where the leading logs (LL) and the next-to-leading logs (NLL), α m s (µ) ln n τ (m ≤ n) in the variable τ = ln(1 + cos φ)/2, have to be resummed. For the fixed-order perturbative calculations, we have used the NLOJET++ package [43,44] and have tested it against the distributions obtained by Madgraph [50]. To achieve numerical stability, we have generated 10 9 DIS events at HERA ( √ s = 314 GeV), allowing us to reach an statistical accuracy of a few percent over most of the phase space.
Being weighted by the product of transverse energies of jets, both the TEEC and ATEEC are expected to be insensitive to the parton distribution functions (PDFs). To quantify this, we use two PDF sets of relatively recent vintage, the CT18 [45], and MMHT14 [46]. The main theoretical uncertainty in the jet physics comes from the scale-dependence, of these the so-called factorization scale µ F enters through the PDFs, and the partonic matrix elements depend essentially on the renormalization scale µ R . Detailed studies done for the inclusive jet and dijet data at HERA show that the µ F -dependence of the cross sections is small, but the µ R -dependence is substantial in the NLO accuracy [47,48]. We study these dependencies in TEEC and ATEEC, following the choice of the nominal scale, µ 0 = E T 2 + Q 2 , as advocated in these papers. The E T donate the average of E T . Varying the scales in the range µ 2 F = (0.5, 2)µ 2 0 , we find that the µ F -dependence is small in the TEEC, not exceeding (1 − 2)% over the cos φ range, but the µ R -dependence is found to be significant. Thus, NNLO improvements are needed to reduce the µ Runcertainty. However, fitting the HERA data on TEEC may also effectively reduce the allowed µ R -range. Finally, we show the sensitivity of the TEEC and ATEEC on the strong coupling constant α s (M Z ), for three representative values α s (M Z ) = 0.108, 0.118, 0.128. With the nominal choice of the scales µ F = µ R = µ 0 , and the current central value of α s (M Z ) = 0.118 [49], we show that the differential distributions TEEC(cos φ) and ATEEC(cos φ) are remarakbly stable perturbatively in both the Q 2 -ranges. This remains to be tested in the NNLO accuracy. Our study presented here makes a good case for using the TEEC in DIS-data as a precision test of perturbative QCD, following similar anayses done for the high energy pp data at the LHC.
The rest of this paper is organized as follows. Section II collects the definitions of TEEC and its asymmetry. Experimental cuts to calculate these functions are stated in this section together with the jet algorithm used and the jet definitions. In Sec. III, we present the numerical results calculated at next-to-leading order in α s and estimate the uncertainty in the shape variables TEEC and ATEEC arising from the different PDFs, and the scale-dependence by varying the scale µ F and µ R . Of these, the µ R -dependence is substantial. Fixing the scale µ R to the nominal value µ 0 , which provides a good fit of the inclusive-jet and dijet data at HERA [47,48], we show the sensitvity of the TEEC and ATEEC on α s (M Z ). A comparison of the LO and NLO results is also presented here. We summarise our results in the last section. A check of the NLOJET++ calculation is shown in Appendix-A at the LO, by using the package MadGraph5 aMC@NLO [50] with the MMHT14 PDF set.

II. TRANSVERSE ENERGY-ENERGY CORRELATION AND ITS ASYMMETRY
In the Breit frame, the transverse energy-energy correlation in γ(q) + p → a + b + X involving hadrons or jets is expressed as: where E T,a and E T,b are transverse energies of two jets or hadrons. The δ-function assures that these hadrons or jets are separated by the azimuthal angle φ, and the cross section σ ′ and Σ ′ indicate kinematic cuts on the integrals, defined later. The second expression is valid for a sample of N hard-scattering multi-jet events, labelled by the index A. The associated asymmetry (ATEEC) is then defined as the asymmetry between the forward (cos φ > 0) and backward (cos φ < 0) parts of the TEEC: Due to the factorization of the amplitudes in QCD, the denominator of the first equation in Eq. (1) dσ γp→a+b+X /dE T can be written as a convolution of the parton distribution functions(PDFs) f q/p (x 1 ), where x 1 is the fractional energy of the proton carried by the parton q, and the parton level cross section. In the leading order, this is given by σ γq→b1b2 . As for numerator, it can also be expressed as the convolution of PDFs with 2 → 3 parton level subprocess, in the leading order, such as γq → qgg. Thus, TEEC is calculated from the following expression: where the symbol ⋆ stands for the convolution. Which processes are included in the calculations of the TEEC depend on the theoretical accuracy. In NLO, this involves 2 → 2, 2 → 3 and 2 → 4 partonic subprocesses. Some representative Feynman diagrams of the subprocess are shown in Fig. 1. In the upper row of Fig  As defined in Eq. (3), the TEEC correlation 1 σ ′ dΣ ′ d cos φ is a normalized variable. In particular, the dependence of the TEEC on the PDFs is compensated to a large extent. Thus, to a good approximation, a factorized result is expected, which can be perturbatively improved by including higher orders. We calculate the TEEC and ATEEC close to experimental conditions used by the HERA experiments H1 and ZEUS, which assume a certain selection criteria based on physical cuts on the kinematic variables. They are defined as follows: The basic DIS kinematic variables x and y = Q 2 /(sx) satisfy Besides, we restrict the range of the pseudo-rapidity in the laboratory frame (η lab ) as The pseudorapidity is related to the polar angle θ, defined with respect to the proton beam direction, by η lab = − ln tan(θ/2). We also use the right-handed co-ordinate system of the H1 collaboration, in which the positive z-axis is in the direction of the proton beam, and the nominal interaction point is located at z = 0.
We calculate the TEEC and ATEEC in the Breit frame used by experiments at HERA. In this frame, transverse energy E T is Lorentz invariant and γp → jjjX is the nontrivial process at the leading order. The cuts for the transverse energy of dijet and trijet events, defined in the Breit frame, are as follows: where the E T 2 and E T 3 denote 1 , respectively. These cuts are consistent with the measurement at HERA [47]. Following the practice in the HERA experimental analysis, we use the k T jetalgorithm [51], where the distance measure of partons (i, j) is given by Here B represents the "beam jet" of the proton: particles with small momenta transverse to the beam axis, and R is the cone-size parameter of the jet which we set to R = 1.0 in our calculation. We use two different PDF sets, CT18 [45] and MMHT14 [46], and explore the uncertainty on the TEEC (cos φ) and ATEEC (cos φ) distributions from these two sets in the next section.
It has become customary to determine the QCD coupling constant at the scale µ = M Z [49]. To determine α s (M Z ) from TEEC (cos φ) and ATEEC (cos φ) in DIS, the cross section can be expressed as: where k denotes a parton (quark or gluon), f k (x, µ F ) is the parton density, and σ k (x, µ F , µ R ) is the partonic cross section, which depends on the renormalization scale µ R and the fatorization scale µ F . The partonic cross section is calculated in perturbative QCD as an expansion in α s : As the µ F -dependence is very mild on TEEC, as shown later, the dominant scale-dependence of the cross section enters through the scale µ R , i.e., from α s (µ R ), which we relate to α s (M Z ) on an event-by-event basis in our simulations. The µ R dependence of α s is given by renormalization group equation In the NLO-calculation, the two-loop β-function is used for transcribing α s (µ) to α s (M Z ) with a certain scale µ which is revelent for the jets defined above. The coupling constant α s (µ) is given as Here, n f is the number of quark flavors, which is determined by the scale µ, we have set n f = 5, and Λ is the QCD parameter, which is determined by the value of α s (M Z ). In the LO-calculation, we set b 1 = 0 in the above expression. In our numerical results, we present TEEC(cos φ) and ATEEC(cos φ) calculated in the LO and NLO for the same value of α s (M Z ), which implies a different value of Λ in the LO and NLO. As already stated, the α s (µ)-dependence enters essentially through µ = µ R . Since µ R is not determined uiquely, there will remain a residual scale-dependence in the differential distributions for TEEC (cos φ) and ATEEC (cos φ). In the next section, we show the dependence of TEEC (cos φ) and ATEEC (cos φ) at HERA on the scales µ F , µ R and α s (M Z ).
Before ending this section, we remark that very recently another shape variable involving the azimuthal angle correlation of the lepton and hadron in DIS process has been proposed and calculated in [52], which is defined as where the sum runs over all hadrons and cos φ ℓa is the cosine of the azimuthal angle between the lepton and the hadron. As seen in the second of the above equation, transverse energy of the lepton drops out of this variable. As opposed to the shape variable TEEC, defined here in Eq. (1) for DIS, as well as the EEC/TEEC variables defined earlier in e + e − annihilation [1,2] and pp collisions [28], which involve (transverse) energy weighted azimutal angle correlations between two jets or hadrons, the shape variable defined in [52] is the azimuthal angle correlation between the lepton and a hadron (or a jet) weighted by the transverse energy of a single hadron (or jet). We emphasize that ℓHT EC(cos φ), defined in [52] and Eq. (13), while interesting in its own right, is a different variable from TEEC.
Lepton-jet correlation in DIS has also been studied in [53], and revisited very recently in [54], where a detailed derivation of the formalism used and a phenomenological study relevant for the jet production at HERA are carried out.

III. RESULTS FOR TEEC (cos φ) AND ITS ASYMMETRY ATEEC (cos φ) IN DIS PROCESS AT HERA
For the numerical results presented here in the LO and NLO accuracy, we have used the program NLOJET++ [43,44]. As a cross check on our calculations, we have also used the program Madgraph to calculate the leading order TEEC and ATEEC functions. To compare with the results obtained using NLOJET++, parton-level events are generated in MadGraph5 aMC@NLO [50] with the MMHT14 PDF set. The distributions obtained from the two packages agree well in both the low-Q 2 (5.5 GeV 2 < Q 2 < 80 GeV 2 ) and high-Q 2 (150 GeV 2 < Q 2 < 1000 GeV 2 ) ranges. The details are given in Appendix A. From now on, we shall work only with the NLOJET++.
We have generated 10 9 events to obtain the NLO results and 10 8 events for the LO results in each of the two Q 2 ranges. This large statistics is required to obtain an accuracy of a few percent, which enables us to meaningfully calculate the various parametric dependences intrinsic to the problem at hand. To check the normalisation, we have calculated the two-jet cross sections at √ s = 314 GeV for the ranges of the DIS variables given in the preceding section and compared them with the corresponding HERA data [47] in Table I. The two-jet events selected for this comparison are defined by the following two bins in < E T 2 and the Q 2 -range given below: Theoretical cross sections are obtained using the CT18 [45] PDFs, the scales set to the values µ R = µ F = E T 2 + Q 2 , and α s (M Z ) = 0.118. The NLO cross sections are in excellent agreement with the HERA data, providing a test of the normalisation.
We restrict cos φ in the range [−0.8, 0.8] to avoid the regions φ ≃ 0 • and φ ≃ 180 • which will involve self-correlations (a = b) and virtual corrections to 2 → 2 processes. In calculating these functions, we use α s (M Z ) = 0.118 and have set the fatorization (µ F ) and the renormalization (µ R ) scales to the following values: µ F = µ R = µ 0 = E T 2 + Q 2 . This scale-setting is discussed in the analysis of the jet-data by the H1 Collaboration [48]. The effect of varying the scale µ F which enters in the PDFs has little effect in the inclusive-and dijet-cross sections [48], which we also find for the TEEC (cos φ) and ATEEC(cos φ), shown later in this section. We quantify the uncertainty on the TEEC (cos φ) and ATEEC (cos φ) from the two input PDFs by the following ratios: They are shown in the lower frames of Fig. 2 and Fig. 3, respectively. We note that in the angular range shown, the pdf-related uncertainty ∆[TEEC(cos φ)] pdf and ∆[ATEEC(cos φ)] pdf is mostly below 10%.
Next, we present the fatorization-scale and the renormalization-scale dependence of the TEEC (cos φ) and ATEEC (cos φ), by fixing the other parameters to their nominal values, and use the MMHT14 pdf set. Fixing µ R = µ 0 , we vary µ F in the range µ 2 F = [0.5, 2]µ 2 0 and show the results of µ F -dependence in Fig. 4 and Fig. 5 for the low-Q 2 range 5.5 GeV 2 < Q 2 < 80 GeV 2 and the high-Q 2 range 150 GeV 2 < Q 2 < 1000 GeV 2 , respectively, for both the LO and the NLO accuracy. The µ F -dependence of TEEC (cos φ) is small, typically about 5%, decreasing for the high-Q 2 range. It is smaller for the asymmetry ATEEC (cos φ), except for the last bin, where it has large statistical error. The The µ R -dependence in the corresponding Q 2 -ranges are shown in Figs. 6 and 7, respectively. Here, we fixed µ F = µ 0 , and varied µ R in the range µ 2 R = [0.5, 2]µ 2 0 . One notices marked improvement in the µ R -dependence from the LO to NLO. However, It is larger than the corresponding µ F -dependence. The µ R -uncertainty on TEEC (cos φ) and ATEEC (cos φ) are plotted in the lower frames of these figures (Figs. 6 and 7) in terms of the ratios ∆[TEEC(cos φ)] µR , and ∆[ATEEC(cos φ)] µR , defined in Eq. (16). Based on these numerical results, we find that the combined uncertaity due to the PDFs, and the µ F and µ R -scales, is at about 10% in the TEEC (cos φ), and smaller in ATEEC (cos φ).
Further reduction in the scale uncertainty requires additional input, which we anticipate from the NNLO improvements as well as from the fits of the HERA data. This is suggested by the detailed NLO-and NNLO-studies done for the inclusive-jet and dijet data at HERA [48], which can be summarized as follows: The effect of varying µ F in the range 10 to 90 GeV on the jet cross sections is small, and this scale can be fixed to a value within this range without risking a perceptible change elsewhere, which is essentially in line what we find in our analyis. The effect of varying the scale µ R is found more significant in the HERA jet-analysis. However, the choice µ R = E T 2 + Q 2 yields a good fit of the jet data in both the NLO and NNLO accuracy. The reduced µ R -dependence in the NNLO accuracy leads to a factor 2 improvement in the accuracy of α s (M Z ). Following [48], we shall fix the scale µ R to its nominal value in studying the sensitivity of TEEC (cos φ) and ATEEC (cos φ) on α s (M Z ).     [49]. They are given in Fig. 8    in the leading order (upper frames) and the next to leading order (lower frames) as in Fig. 4, but for the high-Q 2 range 150 GeV 2 < Q 2 < 1000 GeV 2 at HERA. range) and Fig. 9 (high−Q 2 range) at the LO (upper frame) and NLO accuracy (lower frame). To quantify the α s (M Z )-sensitivity, we define the following ratios:    Renormalization scale dependence of the differential distribution 1/σ ′ dΣ ′ /d(cos φ) and its asymmetry 1/σ ′ dΣ ′asym /d(cos φ) in the leading order (upper frames ) and the next to leading order(lower frames ) varying µ 2 R in the range [0.5, 2] × µ 2 0 , where µ0 is the nominal scale defined in the text, calculated with the MMHT14 PDFs for the low-Q 2 : 5.5 GeV 2 < Q 2 < 80 GeV 2 at HERA. The corresponding µR-dependence is also shown in terms of ∆[TEEC(cos φ)]µ R and ∆[ATEEC(cos φ)]µ R , defined in Eq. (16).
They are shown in the bottom frames in Fig. 8 (low−Q 2 range) and Fig. 9 (high−Q 2 range). Concentrating on the NLO results, we see that both ∆[TEEC(cos φ)] αs and ∆[ATEEC(cos φ)] αs show a marked sensitivity to α s (M Z ), with these ratios reaching as high a value as 0.40. Hence, these shape functions at HERA offer competitive avenues to determine α s (M Z ), and we urge our experimental colleagues to undertake a detailed data analysis of these variables   Renormalization scale dependence of the differential distribution 1/σ ′ dΣ ′ /d(cos φ) and its asymmetry 1/σ ′ dΣ ′asym /d(cos φ) in the leading order (upper frames ) and the next to leading order(lower frames ) as in Fig. 6, but for the high-Q 2 range: 150 GeV 2 < Q 2 < 1000 GeV 2 at HERA.
at HERA. A comparison of the LO and the NLO TEEC(cos φ) and its asymmetry ATEEC (cos φ) at HERA ( √ s = 314 GeV) in the high-Q 2 range and the low-Q 2 range are shown in Fig. 10. These results are obtained for the choice µ F = µ R = µ 0 = E T 2 + Q 2 , α s (M Z ) = 0.118, and MMHT14 set of PDFs. They show that theses correlations are remarkably stable against NLO corrections. We conjecture that NNLO corrections are, likewise, small. This remains to be shown and we hope that our work will stimulate working them out.

IV. SUMMARY
In this paper, we have studied for the first time, the transverse energy energy correlations TEEC (cos φ) and its asymmetry ATEEC (cos φ) in deep inelastic scattering at the electron-proton collider HERA at the center of mass energy √ s = 314 GeV, where φ is the angle in the Breit frame between two jets defined using a transverse-momentum (k T ) jet algorithm. We use NLOJET++ to calculate these functions in the LO and the NLO approximations in  9: Dependence of the differential distribution 1/σ ′ dΣ ′ /d(cos φ) (left) and its asymmetry 1/σ ′ dΣ ′asym /d(cos φ) (right) on the QCD coupling constant αs(MZ ) as in Fig. 8, but for the high-Q 2 range at HERA.
QCD for two ranges in the momentum transfer squared Q 2 . In the LO, these results are checked using the package MadGraph5 aMC@NLO [50] with the MMHT14 PDF set. We show the sensitivity of these functions on the PDFs, factorization (µ F ) and renormalization (µ R ) scales, and on α s (M Z ). With the various cuts in the event generation matched with the ones in the measurements by the H1 collaboration at HERA, these studies are useful in the analysis of the HERA data, including the determination of α s (M Z ) from the shape variables. An NNLO calculation for these shape variables is still lacking. This has the consequence that significant renormalization-scale dependence which enters in the partonic cross sections remains. At the present theoretical FIG. 10: A comparison of the LO and the NLO differential distribution 1/σ ′ dΣ ′ /d(cos φ) (left) and its asymmetry 1/σ ′ dΣ ′asym /d(cos φ) (right) at HERA ( √ s = 314 GeV) in the high-Q 2 range (upper frames) and low-Q 2 range (lower frames), with µF = µR = µ0 = ET 2 + Q 2 and αs(MZ ) = 0.118. accuracy followed in this paper, this may compromise the precision on α s (M Z ). Theoretical precision can be improved by including the NNLO contribution, as shown for the dijet and inclusive jet cross sections in in DIS [55][56][57][58]. However, the scale uncertainty could also be reduced by analysing the HERA data for the shape variables by narrowing the allowed range of µ R for which one gets a good quality fit. This is the case in the analysis of the inclusive-jet and dijet HERA data, in which the choice µ R = E T 2 + Q 2 accounts well the H1 measurements, also in the NLO accuracy [48]. For this choice of the µ R scale, we have shown that the event shape TEEC (cos φ) and its asymmetry are very sensitive to the value of α s (M Z ). We hope that our case-study for the TEEC and ATEEC at HERA, carried out at the NLO accuracy, will help fous on the analysis of the data on these shape vaiables with improved theoretical accuracy.