UvA-DARE ( Digital Academic Repository ) Subjet distributions in deep inelastic scattering at HERA

Subjet distributions were measured in neutral current deep inelastic ep scattering with the ZEUS detector at HERA using an integrated luminosity of 81.7 pb−1. Jets were identified using the kT cluster algorithm in the laboratory frame. Subjets were defined as jet-like substructures identified by a reapplication of the cluster algorithm at a smaller value of the resolution parameter ycut. Measurements of subjet distributions for jets with exactly two subjets for ycut = 0.05 are presented as functions of observables sensitive to the pattern of parton radiation and to the colour eAlso working at Max Planck Institute, Munich, Germany. fNow at Institute of Aviation, Warsaw, Poland. gSupported by the research grant no. 1 P03B 04529 (2005–2008). hThis work was supported in part by the Marie Curie Actions Transfer of Knowledge project COCOS (contract MTKD-CT-2004-517186). iNow at University of Bonn, Germany. jNow at DESY group FEB, Hamburg, Germany. kAlso at Moscow State University, Russia. lNow at University of Liverpool, UK. mOn leave of absence at CERN, Geneva, Switzerland. nNow at CERN, Geneva, Switzerland. oNow at Bologna University, Bologna, Italy. pAlso at Institut of Theoretical and Experimental Physics, Moscow, Russia. qAlso at INP, Cracow, Poland. rAlso at FPACS, AGH-UST, Cracow, Poland. sPartially supported by Warsaw University, Poland. tPartly supported by Moscow State University, Russia. uRoyal Society of Edinburgh, Scottish Executive Support Research Fellow. vAlso affiliated with DESY, Germany. wAlso at University of Tokyo, Japan. xNow at Kobe University, Japan. ySupported by DESY, Germany. zPartly supported by Russian Foundation for Basic Research grant no. 05-02-39028-NSFC-a. aaSTFC Advanced Fellow. abNee Korcsak-Gorzo. acThis material was based on work supported by the National Science Foundation, while working at the Foundation. adNow at University of Kansas, Lawrence, USA. aeAlso at Max Planck Institute, Munich, Germany, Alexander von Humboldt Research Award. afNow at KEK, Tsukuba, Japan. agNow at Nagoya University, Japan. ahMember of Department of Radiological Science, Tokyo Metropolitan University, Japan. aiNow at SunMelx Co. Ltd., Tokyo, Japan. ajAlso at Hamburg University, Inst. of Exp. Physics, Alexander von Humboldt Research Award and partially supported by DESY, Hamburg, Germany. coherence between the initial and final states. Perturbative QCD predictions give an adequate description of the data.


Introduction
Jet production in ep collisions provides a wide testing ground of perturbative QCD (pQCD). Measurements of differential cross sections for jet production [1][2][3] have allowed detailed studies of parton dynamics, tests of the proton and photon parton distribution functions (PDFs) as well as precise determinations of the strong coupling constant, α s . Gluon emission from primary quarks was investigated [4,5] by means of the internal structure of jets; these type of studies gave insight into the transition between a parton produced in a hard process and the experimentally observable jet of hadrons. The pattern of parton radiation within a jet is dictated in QCD by the splitting functions. These functions, P ab (z, µ) with a, b = q or g, are interpreted as the probability that a parton of type b, having radiated a parton of type a, is left with a fraction z of the longitudinal momentum of the parent parton and a transverse momentum squared smaller than µ 2 , where µ is the typical hard scale of the process. The splitting functions are calculable as power series in α s . Thus, the characteristics of jet substructure provide direct access to the QCD splitting functions and their dependence on the scale.
The understanding of jet substructure is also important in the context of jet identification in boosted systems, like hadronic top decays [6] or bb final states at LHC [7]. The first example calls for a direct application of jet substructure, the second requires knowledge about jet substructure to distinguish between single-and double-quark induced jets. This paper presents a study of jet substructure in a more controlled hadronic-type environment than that provided by hadron-hadron colliders.
Jet production in neutral current (NC) deep inelastic scattering (DIS) was previously used to study the mean subjet multiplicity [4] and the mean integrated jet shape [5] with values of α s (M Z ) extracted from those measurements. In the present study, the pattern of QCD radiation is investigated by means of the subjet topology, providing a more stringent test of the pQCD calculations.
In this paper, measurements of normalised differential subjet cross sections for those jets which contain two subjets at a given resolution scale are presented. The measurements were done as functions of the ratio between the subjet transverse energy and that of the jet, E sbj T /E jet T , the difference between the subjet pseudorapidity 1 (azimuth) and that of the jet, η sbj − η jet (|φ sbj − φ jet |), and α sbj , the angle, as viewed from the jet centre, between the subjet with higher transverse energy and the proton beam line in the pseudorapidity-azimuth plane (see Fig. 1). The predictions of pQCD at next-to-leading order (NLO) were compared to the data.

Jets and subjets
The analysis of subjets presented in this paper was performed using the laboratory frame. In this frame, the calculations of the subjet distributions can be performed up to O(α 2 s ), i.e. NLO, with jets consisting of up to three partons. The analysis used events with high virtuality of the exchanged boson, Q 2 ; at low values of Q 2 , the sample of events with at least one jet of high E jet T (E jet T ≫ Q 2 ) is dominated by dijet events. In that case, the calculations include jets consisting of up to only two partons and, therefore, correspond to lowest-order predictions of jet substructure.
The k T cluster algorithm [8] was used in the longitudinally invariant inclusive mode [9] to define jets in the hadronic final state. Subjets [10] were resolved within a jet by considering all particles associated with the jet and repeating the application of the k T cluster algorithm until, for every pair of particles i and j the quantity d ij = min(E T,i , E T,j ) 2 · ((η i − η j ) 2 + (φ i − φ j ) 2 ), where E T,i , η i and φ i are the transverse energy, pseudorapidity and azimuth of particle i, respectively, was greater than d cut = y cut · (E jet T ) 2 . All remaining clusters were called subjets.
The subjet multiplicity depends upon the value chosen for the resolution parameter y cut . Subjet distributions were studied for those jets with exactly two subjets at a value of the resolution parameter of y cut = 0.05. This value of y cut was chosen as a compromise between resolution, size of the hadronisation correction factors and statistics. The effect of the parton-to-hadron corrections on the shape of the subjet distributions becomes increasingly larger as y cut decreases. On the other hand, the number of jets with exactly two subjets decreases rapidly as y cut increases.
Subjet distributions were studied as functions of E sbj T /E jet T , η sbj −η jet , |φ sbj −φ jet | and α sbj . One of the goals of this study was to investigate the extent to which pQCD calculations are able to reproduce the observed distributions. In addition, the dependence of the splitting functions P ab (z, µ) on z can be investigated using the E sbj T /E jet T distribution. The splitting functions at leading order (LO) do not depend on µ but acquire a weak dependence due to higher-order corrections. Such a dependence can be investigated by measuring the subjet distributions in different regions of E jet T or Q 2 . The substructure of jets consisting of a quark-gluon pair (the quark-induced process eq → eqg) or a quark-antiquark pair (the gluon-induced process eg → eqq) are predicted to be different (see Section 8.1). Furthermore, the relative contributions of quark-and gluon-induced processes vary with Bjorken x and Q 2 . The predicted difference mentioned above is amenable to experimental investigation by comparing the shape of the subjet distributions in different regions of x and Q 2 .
Colour coherence leads to a suppression of soft-gluon radiation in certain regions of phase space. The effects of colour coherence between the initial and final states have been studied in hadron-hadron collisions [11]. These effects are also expected to appear in lepton-hadron collisions. For the process eq → eqg, colour coherence implies a tendency of the subjet with lower (higher) transverse energy, E sbj T,low (E sbj T,high ), to have η sbj − η jet > 0 (η sbj − η jet < 0). The variable α sbj , defined in close analogy to the variables used to study colour coherence in hadron-hadron collisions [11], reflects directly whether the subjet with the lower transverse energy has a tendency to be emitted towards the proton beam direction.

Experimental set-up
A detailed description of the ZEUS detector can be found elsewhere [12,13]. A brief outline of the components most relevant for this analysis is given below.
Charged particles were tracked in the central tracking detector (CTD) [14], which operated in a magnetic field of 1.43 T provided by a thin superconducting solenoid. The CTD consisted of 72 cylindrical drift-chamber layers, organised in nine superlayers covering the polar-angle region 15 • < θ < 164 • . The transverse-momentum resolution for full-length tracks can be parameterised as σ(p T )/p T = 0.0058p T ⊕ 0.0065 ⊕ 0.0014/p T , with p T in GeV. The tracking system was used to measure the interaction vertex with a typical resolution along (transverse to) the beam direction of 0.4 (0.1) cm and to cross-check the energy scale of the calorimeter.
The high-resolution uranium-scintillator calorimeter (CAL) [15] covered 99.7% of the total solid angle and consisted of three parts: the forward (FCAL), the barrel (BCAL) and the rear (RCAL) calorimeters. Each part was subdivided transversely into towers and longitudinally into one electromagnetic section and either one (in RCAL) or two (in BCAL and FCAL) hadronic sections. The smallest subdivision of the calorimeter was called a cell. Under test-beam conditions, the CAL single-particle relative energy resolutions were σ(E)/E = 0.18/ √ E for electrons and σ(E)/E = 0.35/ √ E for hadrons, with E in GeV.
The luminosity was measured from the rate of the bremsstrahlung process ep → eγp. The resulting small-angle energetic photons were measured by the luminosity monitor [16], a lead-scintillator calorimeter placed in the HERA tunnel at Z = −107 m.

Data selection
The data were collected during the running period 1998-2000, when HERA operated with protons of energy E p = 920 GeV and electrons or positrons 2 of energy E e = 27.5 GeV, and correspond to an integrated luminosity of 81.7 ± 1.9 pb −1 .
Neutral current DIS events were selected offline using criteria similar to those reported previously [5]. The main steps are given below.
A reconstructed event vertex consistent with the nominal interaction position was required and cuts based on tracking information were applied to reduce the contamination from beam-induced and cosmic-ray background. The scattered-electron candidate was identified using the pattern of energy deposits in the CAL [17]. The energy, E ′ e , and polar angle, θ e , of the electron candidate were also determined from the CAL measurements. The double-angle method [18], which uses θ e and an angle γ that corresponds, in the quark-parton model, to the direction of the scattered quark, was used to reconstruct Q 2 . The angle γ was reconstructed using the CAL measurements of the hadronic final state.
Electron candidates were required to have an energy E ′ e > 10 GeV, to ensure a high and well understood electron-finding efficiency and to suppress background from photoproduction. The inelasticity variable, y, as reconstructed using the electron energy and polar angle, was required to be below 0.95; this condition removed events in which fake electron candidates from photoproduction background were found in the FCAL. The requirement 38 < (E − p Z ) < 65 GeV, where E is the total CAL energy and p Z is the Z component of the energy measured in the CAL cells, was applied to remove events with large initial-state radiation and to reduce further the photoproduction background. Remaining cosmic rays and beam-related background were rejected by requiring the total missing transverse momentum, p miss T , to be small compared to the total transverse energy, E tot The kinematic range was restricted to Q 2 > 125 GeV 2 .
The k T cluster algorithm was used in the longitudinally invariant inclusive mode to reconstruct jets in the measured hadronic final state from the energy deposits in the CAL cells.
The jet algorithm was applied after excluding those cells associated with the scatteredelectron candidate. Jet transverse-energy corrections were computed using the method developed in a previous analysis [5]. Events were required to have at least one jet of E jet T > 14 GeV and −1 < η jet < 2.5. The final sample of 128986 events contained 132818 jets, of which 21162 jets had exactly two subjets at y cut = 0.05.

Monte Carlo simulation
Samples of events were generated to determine the response of the detector to jets of hadrons and the correction factors necessary to obtain the hadron-level subjet cross sections. The hadron level is defined as those hadrons with lifetime τ ≥ 10 ps. The generated events were passed through the Geant 3.13-based [19] ZEUS detector-and trigger-simulation programs [13]. They were reconstructed and analysed applying the same program chain as to the data.
Neutral current DIS events including radiative effects were simulated using the Heracles 4.6.1 [20] program with the Djangoh 1.1 [21] interface to the hadronisation programs. Heracles includes corrections for initial-and final-state radiation, vertex and propagator terms, and two-boson exchange. The QCD cascade is simulated using the colour-dipole model (CDM) [22] including the LO QCD diagrams as implemented in Ariadne 4.08 [23] and, alternatively, with the MEPS model of Lepto 6.5 [24]. The CTEQ5D [25] proton PDFs were used for these simulations. Fragmentation into hadrons is performed using the Lund string model [26] as implemented in Jetset [27,28].
The jet search was performed on the Monte Carlo (MC) events using the energy measured in the CAL cells in the same way as for the data. The same jet algorithm was also applied to the final-state particles (hadron level) and to the partons available after the parton shower (parton level) to compute hadronisation correction factors (see Section 6).

QCD calculations
The O(α 2 s ) NLO QCD calculations used to compare with the data are based on the program Disent [29]. The calculations used a generalised version of the subtraction method [30] and were performed in the massless MS renormalisation and factorisation schemes. The number of flavours was set to five; the renormalisation (µ R ) and factorisation (µ F ) scales were set to µ R = µ F = Q; α s was calculated at two loops using Λ (5) MS = 220 MeV which corresponds to α s (M Z ) = 0.118. The ZEUS-S [31] parameterisations of the proton PDFs were used. The results obtained with Disent were cross-checked by using the program Nlojet++ [32].
Since the measurements refer to jets of hadrons, whereas the QCD calculations refer to jets of partons, the predictions were corrected to the hadron level using the MC samples described in Section 5. The multiplicative correction factor, C had , defined as the ratio of the cross section for subjets of hadrons to that of partons, was estimated with the Lepto-MEPS model, since it reproduced the shape of the QCD calculations better.
The normalised cross-section calculations changed typically by less than ±20% upon application of the parton-to-hadron corrections, except at the edges of the distributions, where they changed by up to ±50%. Other effects not accounted for in the calculations, namely QED radiative corrections and Z 0 exchange, were found to be very small for the normalised cross-section calculations and neglected.
The following theoretical uncertainties were considered (as examples of the size of the uncertainties, average values of the effect of each uncertainty on the normalised cross section as functions of E sbj T /E jet T , η sbj −η jet , |φ sbj −φ jet | and α sbj are given in parentheses): • the uncertainty in the modelling of the parton shower was estimated by using different models (see Section 5) to calculate the parton-to-hadron correction factors (5.6%, 13.2%, 7.6%, 5.3%); • the uncertainty on the calculations due to higher-order terms was estimated by varying µ R by a factor of two up and down (0.01%, 0.46%, 0.58%, 0.34%); • the uncertainty on the calculations due to the choice of µ F was estimated by varying µ F by a factor of two up and down (0.05%, 0.43%, 0.11%, 0.12%); • the uncertainty on the calculations due to those on the proton PDFs was estimated by repeating the calculations using 22 additional sets from the ZEUS analysis [31]; this analysis takes into account the statistical and correlated systematic experimental uncertainties of each data set used in the determination of the proton PDFs (0.07%, 0.18%, 0.12%, 0.05%); • the uncertainty on the calculations due to that on α s (M Z ) was estimated by repeating the calculations using two additional sets of proton PDFs, for which different values of α s (M Z ) were assumed in the fits. The difference between the calculations using these various sets was scaled to reflect the uncertainty on the current world average of α s [33] (0.02%, 0.04%, 0.05%, 0.01%).
These uncertainties were added in quadrature and are shown as hatched bands in the figures.

Corrections and systematic uncertainties
The sample of events generated with CDM, after applying the same offline selection as for the data, gives a reasonably good description of the measured distributions of the kinematic, jet and subjet variables; the description provided by the MEPS sample is somewhat poorer. The comparison of the measured subjet distributions and the MC simulations is shown in Fig. 2.
The normalised differential cross sections were obtained from the data using the bin-by-bin correction method, where N data,i is the number of subjets in data in bin i of the subjet variable A, N had MC,i (N det MC,i ) is the number of subjets in MC at hadron (detector) level, L is the integrated luminosity and ∆A i is the bin width. The MC samples of CDM and MEPS were used to compute the acceptance correction factors to the subjet distributions. These correction factors took into account the efficiency of the trigger, the selection criteria and the purity and efficiency of the jet and subjet reconstruction.
The following sources of systematic uncertainty were considered for the measured subjet cross sections (as examples of the size of the uncertainties, average values of the effect of each uncertainty on the normalised cross section as functions of E sbj T /E jet T , η sbj − η jet , |φ sbj − φ jet | and α sbj are given in parentheses): • the deviations in the results obtained by using either CDM or MEPS to correct the data from their average were taken to represent systematic uncertainties due to the modelling of the parton shower (0.5%, 2.9%, 2.6%, 1.3%); • variations in the simulation of the CAL response to low-energy particles (0.3%, 1.6%, 1.2%, 0.6%).
Other uncertainties, such as those arising from the uncertainty in the absolute energy scale of the jets [1,34], the uncertainty in the simulation of the trigger and the uncertainty in the absolute energy scale of the electron candidate [35], were investigated and found to be negligible. The systematic uncertainties were added in quadrature to the statistical uncertainties and are shown as error bars in the figures.

Results
Normalised differential subjet cross sections were measured for Q 2 > 125 GeV 2 for jets with E jet T > 14 GeV and −1 < η jet < 2.5 which have exactly two subjets for y cut = 0.05. The distribution of the fraction of transverse energy, (1/σ)(dσ/d(E sbj T /E jet T )), is presented in Fig. 3a. It contains two entries per jet and is symmetric with respect to E sbj T /E jet T = 0.5 by construction. This distribution has a peak for 0.4 < E sbj T /E jet T < 0.6, which shows that the two subjets tend to have similar transverse energies.
The η sbj − η jet data distribution is shown in Fig. 3b and also has two entries per jet. The measured cross section has a two-peak structure; the dip around η sbj − η jet = 0 is due to the fact that the two subjets are not resolved when they are too close together. Figure 3c presents the measured normalised cross section as a function of |φ sbj − φ jet |. There are two entries per jet in this distribution. The distribution has a peak for 0.2 < |φ sbj − φ jet | < 0.3; the suppression around |φ sbj − φ jet | = 0 also arises from the fact that the two subjets are not resolved when they are too close together.
The data distribution as a function of α sbj (one entry per jet) increases as α sbj increases (see Fig. 3d). This shows that the subjet with higher transverse energy tends to be in the rear direction. This is consistent with the asymmetric peaks observed in the η sbj − η jet distribution (see Fig. 3b). Figure 4 shows the η sbj − η jet distribution for those jets which have two subjets with asymmetric E sbj , separately for the subjet with higher and lower E sbj T . It is to be noted that since the jet axis is reconstructed as the transverse-energy-weighted average of the subjet axes, the subjet with higher E sbj T is constrained to be closer to the jet axis than that of the lower E sbj T subjet. The measured distributions show that the higher (lower) E sbj T subjet tends to be in the rear (forward) direction. All these observations support the expectation of the presence of colour-coherence effects between the initial and final states and, in particular, the tendency of the subjet with lower E sbj T to be emitted predominantly towards the proton beam direction.

Comparison with NLO QCD calculations
Next-to-leading-order QCD calculations are compared to the data in Figs. 3 and 4. The QCD predictions give an adequate description of the data. However, the data points are situated at the upper (lower) edge of the theoretical uncertainty in some regions of the subjet variables such as E sbj Since the calculations are normalised to unity, the uncertainties are correlated among the points; this correlation gives rise to the pulsating pattern exhibited by the theoretical uncertainties.
The calculation of the cross section as a function of E sbj T /E jet T exhibits a peak at 0.4 < E sbj T /E jet T < 0.6, as seen in the data. The calculations for the η sbj − η jet and α sbj distributions predict that the subjet with higher transverse energy tends to be in the rear direction, in agreement with the data. This shows that the mechanism driving the subjet topology in the data is the eq → eqg and eg → eqq subprocesses as implemented in the pQCD calculations.
To gain further insight into the pattern of parton radiation, the predictions for quark-and gluon-induced processes (see Section 2) are compared separately with the data in Fig. 5. The NLO calculations predict that the two-subjet rate is dominated by quark-induced processes: the relative contribution of quark-(gluon-) induced processes is 81% (19%). The shape of the predictions for these two types of processes are different; in quark-induced processes, the two subjets have more similar transverse energies (see Fig. 5a) and are closer to each other (see Fig. 5b and 5c) than in gluon-induced processes. The comparison with the measurements shows that the data are better described by the calculations for jets arising from a qg pair than those coming from a qq pair.

E jet
T , Q 2 and x dependence of the subjet distributions Figures 6 to 9 show the normalised differential subjet cross sections in different regions of E jet T . Even though the mean subjet multiplicity decreases with increasing E jet T [4], the measured normalised differential subjet cross sections have very similar shapes in all E jet T regions for all the observables considered. This means that the subjet topology does not change significantly with E jet T . This is better illustrated in Fig. 10, where the data for all E jet T regions are plotted together. In particular, it is observed that the maximum of each measured normalised cross section in every region of E jet T occurs in the same bin of the distribution. To quantify the E jet T dependence more precisely, Fig. 11 shows the maximum value of the measured normalised cross section for each observable as a function of E jet T together with the NLO predictions. The spread of the measured maximum values of the normalised cross sections is ±(4 − 6)%. For each observable, the scaling behaviour of the normalised differential subjet cross sections is clearly observed and in agreement with the expectation that the splitting functions depend weakly on the energy scale. The NLO QCD calculations are in agreement with the data and support this observation. Figures 12 to 15 show the normalised differential subjet cross sections in different regions of Q 2 . In this case, it is observed that while the shape of the E sbj T /E jet T distribution does not change significantly with Q 2 , some dependence can be seen in the other observables. For example, the dip in the η sbj − η jet distribution is shallower for 125 < Q 2 < 250 GeV 2 than at higher Q 2 and the shape of the α sbj distribution for 125 < Q 2 < 250 GeV 2 is somewhat different than for the other regions (see Fig. 16). These features of the data are reasonably reproduced by the NLO QCD calculations and understood as a combination of two effects: the fraction of gluon-induced events is predicted to be 32% for 125 < Q 2 < 250 GeV 2 and below 14% for higher Q 2 ; the shape of the normalised cross sections as functions of η sbj −η jet and α sbj changes from the region 125 < Q 2 < 250 GeV 2 to 250 < Q 2 < 500 GeV 2 (see Fig. 17) for quark-and gluon-induced events. It is observed that the maximum of each measured normalised cross section in every region of Q 2 occurs in the same bin of the distribution, except for |φ sbj − φ jet | in the highest-Q 2 region. Figure 18 shows the maximum 3 value of the measured normalised cross section for each observable as a function of Q 2 together with the NLO predictions. The spread of the measured maximum values of the normalised cross sections as functions of E sbj T /E jet T and |φ sbj − φ jet | is ±(3 − 4)%. On the other hand, the measured and predicted maximum values for the normalised cross sections as functions of η sbj − η jet and α sbj exhibit a step-like behaviour between the lowest-Q 2 region and the rest.  Figure 23 shows the data for all x regions plotted together. It is observed that the maximum of each measured normalised cross section in every region of x occurs in the same bin of the distribution, except for |φ sbj − φ jet | in the highest x region. Figure 24 shows the maximum 3 value of the measured normalised cross section for each observable as a function of x. The shape of the E sbj T /E jet T measured distribution does not change significantly with x, whereas some dependence is expected (see Fig. 24a). The dependence of the η sbj − η jet and α sbj distributions with x exhibits features similar to those observed in the study of the Q 2 dependence; in particular, the maximum values (see Figs. 24b and 24d) exhibit a monotonic increase as x increases, which is reasonably reproduced by the calculations. As discussed previously, these features are understood as a combination of two effects: a decrease of the predicted fraction of gluon-induced events from 44% for 0.004 < x < 0.009 to 6% for x > 0.093 and the change in shape of the normalised cross sections for quark-and gluon-induced processes as x increases (see Fig. 25).
To investigate the origin of the change in shape of the normalised differential cross sections between the lowest and higher Q 2 and x regions, LO and NLO calculations were compared. The most dramatic change is observed when restricting the kinematic region to 125 < Q 2 < 250 GeV 2 or 0.004 < x < 0.009 (see Fig. 26); the LO calculation of the η sbj − η jet distribution does not exhibit a two-peak structure as seen in the NLO prediction and in the data. In addition, the LO calculation of the α sbj distribution peaks at α sbj ∼ π/2 in contrast with the NLO prediction and the data. This proves that the NLO QCD radiative corrections are responsible for these variations in shape and necessary for describing the data.
In summary, while the shapes of the normalised differential cross sections show only a weak dependence on E jet T , their dependence on Q 2 and x have some prominent features at low Q 2 or x. The weak dependence on E jet T is consistent with the expected scaling behaviour of the splitting functions; however, the restriction to low Q 2 or x values demonstrates that the NLO QCD radiative corrections are important there. The NLO QCD calculations, which include the two competing processes eq → eqg and eg → eqq and radiative corrections, adequately reproduce the measurements.

Summary
Normalised differential subjet cross sections in inclusive-jet NC DIS were measured in ep collisions using 81.7 pb −1 of data collected with the ZEUS detector at HERA. The cross sections refer to jets identified in the laboratory frame with the k T cluster algorithm in the longitudinally invariant inclusive mode and selected with E jet T > 14 GeV and −1 < η jet < 2.5. The measurements were made for those jets which have exactly two subjets for y cut = 0.05 in the kinematic region defined by Q 2 > 125 GeV 2 .
The cross sections were measured as functions of E sbj T /E jet T , η sbj −η jet , |φ sbj −φ jet | and α sbj . The data show that the two subjets tend to have similar transverse energies and that the subjet with higher transverse energy tends to be in the rear direction. This is consistent with the effects of colour coherence between the initial and final states, which predict that soft parton radiation is emitted predominantly towards the proton beam direction.
An adequate description of the data is given by NLO QCD calculations. This means that the pattern of parton radiation as predicted by QCD reproduces the subjet topology in the data. Furthermore, the subjet distributions in the data are better described by the calculations for jets arising from a quark-gluon pair.
The normalised cross sections show a weak dependence on E jet T , in agreement with the expected scaling behaviour of the splitting functions. By restricting the measurements to low Q 2 or x values, significant differences in shape are observed, which can be primarily attributed to NLO QCD radiative corrections.  (c) (d) Figure 11: Maximum of the measured normalised differential (a) E sbj T /E jet T , (b) η sbj − η jet , (c) |φ sbj − φ jet | and (d) α sbj subjet cross sections (dots) for jets with E jet T > 14 GeV and −1 < η jet < 2.5 which have two subjets for y cut = 0.05 in the kinematic region given by Q 2 > 125 GeV 2 as a function of E jet T . For comparison, the NLO predictions for quark-(dotted histograms) and gluon-induced (dot-dashed histograms) processes are also shown separately. Other details are as in the caption to Fig. 3. (c) (d) Figure 17: Predicted normalised differential subjet cross sections (solid histograms) for jets with E jet T > 14 GeV and −1 < η jet < 2.5 which have two subjets for y cut = 0.05 in the kinematic region given by Q 2 > 125 GeV 2 as functions of (a,c) η sbj − η jet and (b,d) α sbj in different regions of Q 2 . The NLO predictions for quark-(dotted histograms) and gluon-induced (dot-dashed histograms) processes are also shown separately. (c) (d) Figure 18: Maximum of the measured normalised differential (a) E sbj  (c) (d) Figure 24: Maximum of the measured normalised differential (a) E sbj T /E jet T , (b) η sbj − η jet , (c) |φ sbj − φ jet | and (d) α sbj subjet cross sections (dots) for jets with E jet T > 14 GeV and −1 < η jet < 2.5 which have two subjets for y cut = 0.05 in the kinematic region given by Q 2 > 125 GeV 2 as a function of x. For comparison, the NLO predictions for quark-(dotted histograms) and gluon-induced (dot-dashed histograms) processes are also shown separately. Other details are as in the caption to Fig. 3. (c) (d) Figure 25: Predicted normalised differential subjet cross sections (solid histograms) for jets with E jet T > 14 GeV and −1 < η jet < 2.5 which have two subjets for y cut = 0.05 in the kinematic region given by Q 2 > 125 GeV 2 as functions of (a,c) η sbj − η jet and (b,d) α sbj in different regions of x. The NLO predictions for quark-(dotted histograms) and gluon-induced (dot-dashed histograms) processes are also shown separately. (c) (d) Figure 26: Measured normalised differential subjet cross sections (dots) for jets with E jet T > 14 GeV and −1 < η jet < 2.5 which have two subjets for y cut = 0.05 in restricted Q 2 and x regions as functions of (a,c) η sbj − η jet and (b,d) α sbj . The NLO (solid histograms) and LO (dashed histograms) calculations are also shown. The hatched bands represent the NLO theoretical uncertainty.