Three- and Four-jet Production at Low x at HERA

Three- and four-jet production is measured in deep-inelastic $ep$ scattering at low $x$ and $Q^2$ with the H1 detector using an integrated luminosity of $44{.}2 {\rm pb}^{-1}$. Several phase space regions are selected for the three-jet analysis in order to study the underlying parton dynamics from global topologies to the more restrictive regions of forward jets close to the proton direction. The measurements of cross sections for events with at least three jets are compared to fixed order QCD predictions of ${\mathcal{O}}(\alpha_{\rm s}^2)$ and ${\mathcal{O}}(\alpha_{\rm s}^3) $ and with Monte Carlo simulation programs where higher order effects are approximated by parton showers. A good overall description is provided by the ${\mathcal{O}}(\alpha_{\rm s}^3) $ calculation. Too few events are predicted at the lowest $x \sim 10^{-4}$, especially for topologies with two forward jets. This hints to large contributions at low $x$ from initial state radiation of gluons close to the proton direction and unordered in transverse momentum. The Monte Carlo program in which gluon radiation is generated by the colour dipole model gives a good description of both the three- and the four-jet data in absolute normalisation and shape.


Introduction
The HERA electron-proton collider has extended significantly the available kinematic range for tests of Quantum Chromodynamics. The high centre of mass energy of 319 GeV allows for deep-inelastic scattering (DIS) at a large negative four momentum transfer squared Q 2 ≥ 4 GeV 2 on partons which carry a very small fraction x of the proton momentum down to values of 10 −4 . This is the domain of high parton densities in the proton dominated by gluons and sea quarks. In addition, DIS at low x corresponds to scattering at high γ * p centre of mass energies and is therefore intimately linked to the high energy behaviour of QCD.
Many of the available calculations for DIS processes make use of collinear factorisation [1]: the cross sections are expressed as a convolution of hard partonic subprocesses with proton parton density functions (PDFs). The latter describe the probabilities to find partons in the proton which carry a fraction x of the proton momentum. The separation of the calculation into two pieces is specified by the factorisation scale µ 2 f : initial state radiations from the proton with virtualities above this scale are treated in the hard partonic part while those below are absorbed in the PDFs. For inclusive DIS, Q 2 provides the natural scale for µ f , i.e. µ 2 f = Q 2 . The evolution of the PDFs with µ 2 f is generally described by the DGLAP [1] equations. To leading logarithmic accuracy this is equivalent to the exchange of a parton cascade, with the exchanged partons strongly ordered in virtuality up to Q 2 . For low x this becomes approximately an ordering in k T , the transverse momentum of the partons in the cascade as shown in figure 1. However, at low x the collinear factorisation scheme may break down. The DGLAP leading logarithmic approximation neglects topologies of gluon radiation with unordered k T . These appear in the full perturbative expansion as (log 1/x) n terms which are naturally expected to become large at low x. At very low values of x it is believed that the theoretically most appropriate description is given by the BFKL evolution equations [2], which resum large logarithms of 1/x. The BFKL resummation imposes no restriction on the ordering of the transverse momenta within the parton cascade. Compared to the DGLAP approximation more gluons with sizable transverse momentum are emitted near the proton direction, referred to in the following as the forward direction. This should lead to a significantly increased rate of forward jets [3]. A promising approach to parton evolution at low and larger values of x is given by the CCFM [4] evolution equation, which, because it uses angular-ordered parton emission, is equivalent to the BFKL ansatz for x → 0, while reproducing the DGLAP equations at large x.
Higher order calculations in the collinear factorisation scheme can also improve the treatment for the low x region, since the log 1/x terms can be treated up to the given order of α s . It is an interesting question how well such fixed order approximations can work or if one would still need a full log 1/x resummation to all orders as provided by the BFKL approach. For inclusive DIS ep → eX, QCD analyses (see e.g. [5][6][7][8][9]) were performed using collinear factorisation calculations in next-to-leading order (NLO) O(α s ) 1 and/or next-to-next-to-leading (NNLO) O(α 2 s ) for the hard subprocess with parton densities matched to that order. These calculations are able to describe the inclusive DIS data from HERA and fixed target experiments over a large range in Q 2 down to Q 2 ∼ 2 GeV 2 from the largest to the lowest covered x values.
Final states with jets in DIS are an ideal tool to investigate low x dynamics: the jets can be used to tag higher order processes and furthermore provide direct access to the outgoing hard partons. The H1 and ZEUS measurements of dijet production at low x [10][11][12] and of inclusive forward jet production [13][14][15] show that the leading order (LO) O(α s ) calculations based on DGLAP greatly underestimate the data. NLO O(α 2 s ) calculations can account for some of the LO deficiencies, but the description remains unsatisfactory at low x and Q 2 . In the present paper events with at least three-or four-jets in the final state are investigated. In contrast to inclusive jets and dijets, three-and four-jet final states require the radiation of at least one and two hard gluons respectively in addition to the qq pair from the dominating hard boson-gluonfusion scattering process γ * g → qq (see figure 2). Therefore three-and four-jet processes are ideally suited to study the gluon emissions and the underlying parton dynamics in the proton. Three-jet cross sections in DIS have been measured previously, both by the H1 [16] and the ZEUS [17,18] collaborations. In these analyses the leading jets were required to have a large transverse momentum of at least about 7 GeV. All measured cross sections were found to be well described by NLO O(α 3 s ) predictions in the collinear factorisation scheme. This paper presents a new measurement of three-jet production. The analysis is performed in an extended phase space, covering jets with low transverse momenta down to 4 GeV, and based on a three times larger luminosity than used in the previous H1 publication [16]. The analysis reaches values of x as low as x = 10 −4 . In addition, cross sections for events with at least four jets are measured for the first time in DIS. The data are compared with the NLOJET++ [19] fixed order calculations in the collinear factorisation scheme. This program provides LO O(α 2 s ) predictions for the three-jet case. In addition NLOJET++ is the only available program which provides perturbative calculations for jet cross sections in hadronic collisions to O(α 3 s ) accuracy. This corresponds to NLO and LO precision for the three-and four-jets cases, respectively.
The two Monte Carlo generators RAPGAP [20] and DJANGOH [21], which were able to describe reasonably well inclusive forward jet and dijet production at low x, are also tested.
The sensitivity to deviations from the DGLAP approach may be increased by selecting kinematic regions where gluon radiation is suppressed for this approximation. This is the case for events with a hard forward jet and a large separation in rapidity to a central parton system. Two different subsets of the inclusive three-jet sample are studied: one sample with one forward jet and two central jets and one with two forward jets and one central jet. Figure 2 shows two examples of DIS processes dominating the production of three or more jets at low x. The diagrams contribute to order α 2 s and α 3 s to the cross section, respectively. The radiated gluons are predominantly emitted in the forward direction whereas the quarks from the hard scattering process are mostly central.  Figure 2: Examples of leading order (left) and next to leading order (right) diagrams for three-jet production in DIS at HERA with one and two radiated gluons, respectively.

Kinematics and Measurement Observables
The kinematic variables which describe the hard electron-quark scattering process are the negative four momentum transfer squared Q 2 = −q 2 = (k − k ′ ) 2 of the exchanged virtual photon (γ * ), the Bjørken variable x = Q 2 /(2pq), and the inelasticity y = (qp)/(kp), where k, k ′ , p and q denote the four momenta of the incoming and outgoing positron, the incoming proton and the exchanged photon, respectively. The three variables are related by Q 2 = xys, where s denotes the fixed ep centre of mass energy squared. Jets are defined in this analysis in the γ * p centre of mass system. The observables used to characterise the jets are their transverse momentum p * T in the γ * p centre of mass frame and their pseudorapidity η in the laboratory system. The topology of a three-jet system is fully specified by the following four canonical variables [23]: the scaled energy of the jets ) and the two three-jet angles θ ′ and ψ ′ as defined in figure 3. These variables are measured in the three-jet centre of mass frame. Four-jet events have additional degrees of freedom. The two jets with the lowest dijet mass are combined in order to use the same variables as in the three-jet case.
Three−jet Rest Frame 1 3 2 Figure 3: Definition of the angles θ ′ and ψ ′ in the three-jet rest system [23]. The 3-vector p ′ beam is defined by p ′ beam = p ′ in, 1 − p ′ in, 2 where p ′ in, 1 and p ′ in, 2 are the 3-momenta of the two incoming interacting particles in the three-jet centre of mass frame. They are sorted with respect to their energy in the laboratory frame: E in, 1 > E in, 2 . The incoming interacting particles are the exchanged virtual boson and the parton from the proton side (predominantly a gluon). The latter is assumed to move parallel with the proton and to carry a fraction of its momentum reconstructed as x gluon = x · (1 + (ŝ/Q 2 )), whereŝ denotes the squared centre of mass energy of the three-jet system.

QCD Predictions
The RAPGAP [20] and DJANGOH [21] Monte Carlo event generator programs are used in this analysis to estimate the corrections that must be applied to the data for the finite acceptance, efficiency and resolution of the detector. The two programs are also used to provide predictions that can be compared with the data. Both programs generate hard QCD 2→2 subprocesses (e.g. γ * g → qq) which are convoluted with the CTEQ5L [24] set of parton distributions for the proton. The factorisation and renormalisation scales are set to µ 2 f = µ 2 r = Q 2 for DJANGOH and µ 2 f = µ 2 r = Q 2 +p 2 T for RAPGAP, wherep T is the transverse momentum of the outgoing hard partons. RAPGAP includes resolved photon processes using the SaS 2D photon parton distribution functions [25], which were found to give a good description of the effective photon structure function as measured by H1 [26]. Higher order QCD effects that produce further hard outgoing partons are generated in both RAPGAP and DJANGOH by parton showers: in RAPGAP the showers are ordered in the transverse momenta (k T ) of the emissions, according to the DGLAP leading log Q 2 approximation. DJANGOH uses the colour dipole model (CDM) [22], in which partons are generated by colour dipoles, spanned between the partons in the cascade. Since the dipoles radiate independently, there is no k T ordering. For the hadronisation, the Lund string fragmentation [27] is used both for RAPGAP and DJANGOH. QED radiative corrections are applied in DJANGOH using the HERACLES [28] program and are neglected for RAPGAP. The DJANGOH predictions are referred to as CDM in the following.
Fixed order QCD predictions at parton level are calculated using the NLOJET++ [19] program, which is able to predict three-jet cross sections in LO O(α 2 s ) or NLO O(α 3 s ) and four-jet cross sections in LO O(α 3 s ). The renormalisation and factorisation scales are set to with m = 3 for the three-jet and m = 4 for the four-jet cross sections and N jet being the number of jets fulfilling the applied jet cuts. The value of α s (M Z ) is fixed to 0.118 and the CTEQ6M [7] proton parton density parameterisations are used. The NLOJET++ parton level cross sections are corrected bin-by-bin for hadronisation effects using the CDM simulation as discussed in detail in section 4.3. Two uncertainties are considered for the NLOJET++ cross sections: The uncertainty due to missing higher orders is estimated by recalculating the cross sections with the scales µ f and µ r varied by a common factor of 2 or 0.5. Hadronisation uncertainties are estimated by determining the corrections to the hadron level alternatively with RAPGAP and taking 50% of the difference between the corrections from CDM and RAPGAP as systematic error.

The H1 Detector
A detailed description of the H1 detector can be found in [29]. Here, a brief account of the components most relevant to the present analysis is given. The H1 coordinate system convention defines the outgoing proton beam direction as the positive z axis, also referred to as the 'forward' direction. The polar angle θ is defined with respect to this direction. The pseudorapidity is given by η = − ln tan(θ/2).
The central ep interaction region is surrounded by two large concentric drift chambers (CJCs), operated inside a 1.16 T solenoidal magnetic field. Charged particles are measured in the pseudorapidity range −1.5 < η < 1.5 with a transverse momentum resolution of σ(p T )/p T ≃ 0.005 p T / GeV ⊕ 0.015. Two additional drift chambers (CIZ, COZ) complement the CJCs by precisely measuring the z coordinates of track segments and hence improve the determination of the polar angle. Multi-wire proportional chambers (MWPC) provide fast signals for triggering purposes.
The data sample of this analysis was collected using a trigger which requires the scattered positron to be measured in the SpaCal, at least one high transverse momentum track (p T > 800 MeV) to be reconstructed in the central tracking chambers and an event vertex to be identified by the MWPCs. The trigger efficiency is higher than 85% for the whole analysis phase space.
The ep luminosity is measured via the Bethe-Heitler Bremsstrahlung process ep → ep γ, the final state photon being detected in a crystal calorimeter at z = −103 m.

Event Selection and Kinematic Reconstruction
A detailed account of this analysis can be found in [32]. The data used in this analysis were taken in the 1999 and 2000 running periods, in which HERA collided protons with an energy of 920 GeV with 27.5 GeV positrons, corresponding to a centre of mass energy of √ s = 319 GeV.
The integrated luminosity of the data is 44.2 pb −1 . DIS events are preselected requiring a scattered positron measured in the SpaCal with an energy E e > 9 GeV. The polar angle θ e of the scattered positron is determined from the cluster position in the SpaCal and the z position of the event vertex reconstructed with the central tracking chambers. The observables y, Q 2 and x are derived from the electron kinematics where E e, beam is the positron beam energy. The kinematic range is chosen to be 5 GeV 2 < Q 2 < 80 GeV 2 , 0.1 < y < 0.7, 10 −4 < x < 10 −2 and 156 • < θ e < 175 • .
The hadronic system, containing the jets, is measured with the LAr and SpaCal calorimeters and the central tracking system. Calorimeter cluster energies and track momenta are combined using algorithms which avoid double counting [33]. Jets are formed from the hadronic final state particles boosted to the γ * p rest frame. The inclusive k T cluster algorithm [34] is applied with a separation parameter of 1.0. The p T weighted recombination scheme is used in which the jets are treated massless. The jets are ordered with respect to their transverse momentum in the γ * p rest frame (p * T i > p * T i+1 ). Only jets with a transverse momentum p * T i of at least 4 GeV and a pseudorapidity in the range −1 < η i < 2.5 are considered for the analysis. The latter cut ensures the jets to lie well within the acceptance of the LAr calorimeter. At least three jets are required which fulfil these cuts. It is demanded in addition that p * T 1 + p * T 2 > 9 GeV. The applied cuts ensure a good correlation between jets at detector level and hadron or parton level and allow for comparison of the data to the NLO O(α 3 s ) calculation. In addition to the above selections one of the three leading p * T jets has to lie in the central region −1 < η i < 1.3. This ensures a good trigger efficiency. After all cuts, 38400 events are selected with at least three jets. 5900 of these events contain at least four jets.

Cross Section Measurement
The kinematic region for which the cross sections are measured is given in table 1. All cross sections are given as bin-averaged differential cross-sections defined at the level of stable hadrons. Therefore the data are corrected for all detector effects, using Monte Carlo simulations. For each generated event the response of the H1 detector is simulated in detail including trigger effects. The events are then subjected to the same reconstruction and analysis programs as the data. For each measurement bin a correction factor is calculated as the ratio of simulation entries at stable hadron level to that at detector level. The same inclusive k T algorithm is applied at the hadron and detector levels. The detector correction factors are determined using the CDM simulation which is found to give a better description of the jet topologies than RAPGAP. The Monte Carlo events are weighted in a few variables to adjust their kinematic distributions to the data. These variables are the p * T of the leading jet, η 1 − η 2 , η 1 + η 2 and Q 2 . After weighting the simulations provide a reasonable description of the shapes of all data distributions. The detector correction factors have been studied in detail for all distributions. They vary between 0.6 and 1.2 for events with at least three-jets (0.4 and 1.2 for events with at least four-jets) and show a smooth behaviour. Further small corrections are applied to the data to take QED radiative effects into account. The data are corrected to the QED Born level using the CDM simulation. A correction factor is determined for each measurement bin separately.
For comparison with the data, the fixed order NLOJET++ parton level calculations are corrected to the stable hadron level by application of hadronisation correction factors c had . These corrections are estimated bin-by-bin using the weighted CDM simulation. Jets are obtained at the parton level using the inclusive k T algorithm, both in NLOJET++ and CDM. For CDM the algorithm is applied to the partons after the parton showering step. As just mentioned, the detector and hadronisation corrections are calculated using the weighted Monte Carlo simulation events. However, the unweighted Monte Carlo predictions are compared to the data, as will be shown in section 5.
The correlations between the jets at the different levels have been studied in detail using Monte Carlo simulated events. According to CDM, for the phase space given in table 1, 73-85% of the selected detector level jets can be associated with a hadron level jet within a cone ∆R = (∆η 2 + ∆φ 2 ) ≤ 0.4 around the detector level jet and with a hadron level transverse momentum of p * T ≥ 1.5 GeV. This fraction of "matched" detector jets decreases to 65-75% at the parton level. The lowest matching fractions are observed for the more forward jet pseudorapidities. These migrations dilute the interpretability of the data in terms of the underlying partonic dynamics and must be well controlled. The fraction of unmatched jets observed in RAPGAP agrees with CDM to better than 30% everywhere. Taking the differences from RAP-GAP and CDM as the absolute uncertainty of the number of unmatched jets and assuming this number to directly propagate into the measured cross sections a maximal possible cross section error of 25% is derived. This possible error will be only considered in the discussion of the results (section 5), whenever a large excess of data over NLO prediction is observed. It is not included in the standard determination of systematic cross section errors which is described in the following.

Estimate of Systematic Errors
The errors of the measured differential cross sections are separated into statistical errors of the data δ stat , additional uncorrelated errors δ uncorr , accounting for the statistical errors of the Monte Carlo samples used to determine the various correction factors, and systematic errors. The latter are separated into two contributions: a global normalisation error δ norm and a correlated error δ corr which affects only the shape of the cross section distributions.
The effects of systematic uncertainties on the cross sections are evaluated by applying appropriate variations to the Monte Carlo simulations. The following sources of error are considered: • LAr hadronic energy scale: The absolute hadronic energy scale of the LAr calorimeter is known to 4% accuracy. This is the dominating uncertainty for the determination of the energy of the jets studied in this analysis.
• SpaCal electromagnetic energy scale: The energy of the scattered positron is known within a 2% uncertainty.
• Positron angle measurement: The uncertainty in the measurement of the polar angle of the scattered positron is 1 mrad.
• Track contribution to combined objects: The uncertainty of this contribution is estimated by varying the momenta of all contributing tracks by ±3%.
• Trigger efficiency: The simulated trigger efficiencies are compared with the efficiencies determined from data, using monitor trigger samples. Agreement is found within 3%.
• Luminosity measurement: The measurement of the integrated luminosity is accurate within 1.5%.
• Radiative correction: The uncertainty of the radiative correction factors is estimated to be 2% [35].  Table 2: Sources of correlated systematic uncertainties and the resulting errors on the cross sections for events with at least three jets. The first column contains the error source and the second the range in which the quantity was varied to account for its uncertainty. The two remaining columns give the typical correlated (δ corr ) and global normalisation (δ norm ) uncertainties on the cross sections which arise from each error source.
• Model uncertainty: The cross sections are corrected to hadron level using the weighted CDM simulation. The uncertainty of these corrections is estimated by calculating the correction factors with the weighted RAPGAP simulation and taking 50% of the difference to CDM as systematic uncertainty.
Typical values for the correlated uncertainties and the global normalisation error on the cross sections for events with at least three jets are given in table 2. The systematic errors are dominated by the LAr hadronic energy scale. The second largest contribution stems from the model uncertainty. The total global normalisation error is +16 −14 %. The systematic uncertainties for the cross sections of events with at least four jets are found to be of similar or somewhat larger size as those for events with at least three-jets; the total global normalisation error is +22 −19 %.

Results
The measured cross sections are shown in the figures 4 to 12 and listed in the tables 3 to10.

Cross Sections for Events with at Least Three Jets
Differential cross sections are presented in figures 4 to 6 for events with three or more jets as a function of the number of jets (N jet ), the Bjørken variable x, the pseudorapidities of the three jets and the variables characterising the topology of the three jets in the three-jet centre of mass frame (scaled jet energies and three-jet angles). The kinematic range for which the cross sections are determined is specified in table 1. The figures also show the predictions of the NLOJET++ fixed order QCD prediction in LO O(α 2 s ) and NLO O(α 3 s ), corrected for hadronisation effects. The theory error, including scale variations and hadronisation correction uncertainties added in quadrature, exceed the measurement uncertainty. Figure 4 (left) shows the distribution of the number of jets found in the selected events which extends up to N jet = 6. For this distribution the predictions of the NLO O(α 3 s ) calculation and of the two Monte Carlo programs RAPGAP and CDM are shown. The NLO O(α 3 s ) prediction agrees for N jet = 3 and underestimates the rate of events with 4 jets by a factor 2.6. It does not produce any events with more than four jets as expected.
The differential cross sections shown in figures 4 to 6 are not described by the LO O(α 2 s ) QCD predictions neither in shape nor in magnitude. The main discrepancies are seen at low x and for forward jets (large positive η) where far too few events are predicted. The NLO O(α 3 s ) prediction improves the situation considerably in all regions where deviations from LO are observed. A similar improvement was already noticed in the previous H1 three-jet analysis [16], which was restricted to the phase space of large invariant three-jet masses above 25 GeV. For that phase space the O(α 3 s ) calculation could still describe the three-jet data down to the smallest x = 10 −4 . However, for the present analysis without such mass cuts applied, at very small x < 2 · 10 −4 the calculation undershoots the data, which lies approximately at the upper edge of the total theoretical error band, by a factor of about 0.6. In the complementary region x > 2 · 10 −4 , the description is reasonable; this was also observed in the recent ZEUS multijet analysis [18] which was restricted to this phase space. In summary, a large deficit of the NLO prediction persists only at low x and for forward jets. This is the kinematic region where unordered gluon radiation is expected to enhance the jet production [3]. The shapes of the cross section distributions for the three-jet topological variables shown in figure 6 are all well described by the NLO O(α 3 s ) prediction, only the global normalisation of the calculation is somewhat too low.

Forward Jet Subsample
The observed excess of data versus QCD predictions in the region of forward jet rapidities and at low x is investigated here in further detail. The sample of selected events with at least three jets fulfilling the criteria presented in table 1 is reduced to a subsample by requiring that at least one of the three leading jets be forward and carry a large proton momentum fraction: Here E jet is the jet energy in the laboratory frame and E p, beam = 920 GeV the proton beam energy. Further requirements are applied to obtain two separate subsamples. In the sample with one forward jet and two central jets (f+2c) the other two leading jets are required to lie in the central pseudorapidity range −1 < η jet < 1. In the sample with two forward jets and one central jet (2f+c) it is demanded that one of the three leading jets is a central jet with −1 < η jet < 1 and the remaining leading jet must be in the more forward region η jet > 1 (for this second forward jet no cut is applied on x jet ). The fraction of jets due to gluon radiation is expected to be larger for forward jets than for central jets. This is confirmed by a study of the parton composition of three-jets in the CDM simulation. Therefore the f+2c sample will have many events with only a single radiated gluon (as for the left diagram in figure 2) while the 2f+c selection will have a larger fraction of events with two radiated gluons (as for the right diagram in figure 2). Cross section measurements as a function of x, η 1 and p * T 1 for the f+2c and 2f+c samples are presented in figures 7 and 8. The additional global normalisation errors of +19 −14 % for the f+2c sample and +18 −15 % for the 2f+c selection are not shown in the figures. The fixed order NLO O(α 3 s ) prediction gives a rather good description of the f+2c sample. The step from LO O(α 2 s ) to NLO O(α 3 s ) improves the agreement mainly at very low x < 2 · 10 −4 , where a remaining deficiency of ∼30% is observed. Towards larger x values both the LO and the NLO calculation fall off somewhat less steeply than the data and are too high for the largest covered x values x > 2 · 10 −3 . For the 2f+c selection an even more dramatic change is observed at low x from O(α 2 s ) to O(α 3 s ): the discrepancy at x < 2 · 10 −4 is reduced from a factor of 10 to 2.6. The large remaining deficiency exceeds the combined error of prediction and data and is thus highly significant. It can also not be explained by a possible additional maximal cross section error of 25% due to detector jets which cannot be matched with hard partons (as discussed in the last paragraph of section 4.3). This data excess provides a strong hint for missing higher order QCD corrections, i.e. beyond O(α 3 s ), in this forward gluon radiation dominated phase space. Note that for processes with two radiated gluons, the O(α 3 s ) calculation can only provide a leading order perturbative estimate.
Excesses which are probably related to the one reported here were observed in the forward jet analyses from H1 [13] and ZEUS [15]. In these analyses the topologies of three jets were investigated for events containing a dijet system in addition to a forward jet. The O(α 3 s ) predictions were found to undershoot the data in the region where all three jets tend to go forward. However, in these analyses the data were either integrated over a larger x range or restricted to somewhat larger x values, which might explain why the excesses are less prominent than observed in the present measurement.

Monte Carlo Program Predictions for Events with at Least Three Jets
The cross sections for events with at least three jets are compared to predictions from RAPGAP and CDM. The jet multiplicity shown in figure 4 (left) is described well by CDM while RAP-GAP falls off too steeply. The overall normalisation of the RAPGAP and CDM predictions is found to be too low by 55% and 5% respectively. In the following the Monte Carlo predictions are normalised to the total measured cross section in order to compare only the shapes of the cross sections. Figures 9 and 10 show the comparison as a function of the Bjørken scaling variable x, the difference of the pseudorapidity of the two leading p * T jets (∆η = η 1 − η 2 ), the variables p * T 1 , X ′ 2 and the two three-jet angles cos θ ′ and cos ψ ′ . RAPGAP fails to describe the x and ∆η distributions. CDM on the other hand gives a very good description of almost all observables besides p * T 1 , where it predicts too many high momentum jets (p * T 1 > 15 GeV). The three-jet angular distribution cos ψ ′ is also described rather poorly by both RAPGAP and CDM. A separate check of the cross sections for p * T 1 > 20 GeV reveals significant deviations to CDM for the shapes of various distributions, especially for x and η 1 , as presented in figure 11. In this domain RAPGAP describes these distributions well.

Cross Sections for Events with at Least Four Jets
A subsample of events with four or more jets is also studied. All selection criteria of the threejet sample listed in table 1 have to be fulfilled. In addition, at least one more jet has to be found which satisfies the standard jet cuts p * T > 4 GeV and −1 < η < 2.5. As already mentioned the NLOJET++ O(α 3 s ) calculation can only provide a LO prediction for the final state with four jets which is by far too low as can be seen in figure 4 (left). Thus in the following the comparisons of the measured four jet cross sections are restricted to the CDM and RAPGAP predictions, where parton showers approximate higher orders and can lead to large jet multiplicities. The total cross section predicted by CDM for events with four or more jets agrees well with the data while RAPGAP is too low by a factor of ∼ 2.9, as shown in figure 4 (left). Differential cross sections as a function of p * T 1 , η 1 − η 4 , X ′ 2 and cos θ ′ for events with at least four jets are shown in figure 12 and compared to the predictions by the two Monte Carlo generators normalised to the data. RAPGAP fails to describe the shapes of the differential distributions, again with the exception of the momentum distributions of the jets. CDM on the other hand disagrees with the data in the p * T distributions but describes the scaled energies of the four jets correctly. It also gives a very good description of all other distributions.

Summary
This paper presents a new measurement of three-jet production in DIS at low x and Q 2 . The measurement is carried out in an extended kinematic phase space covering lower jet transverse momenta compared to previous three-jet analyses. Very small x values are reached down to x = 10 −4 . The first measurement of four-jet production in DIS is also presented. Three-and four-jet final states require the radiation of at least one respectively two hard gluons from the initial state proton, in addition to the qq pair from the dominating hard boson-gluon-fusion scattering process γ * g → qq, and are therefore well suited to study parton dynamics at small x.
The measurements are compared with the NLOJET++ [19] fixed order QCD calculations. A remarkable result of the present analysis is the success of the next-to-leading order O(α 3 s ) calculation for the cross sections of events with at least three jets. The inclusion of diagrams with two radiated gluons improves dramatically the agreement with the data compared to the O(α 2 s ) prediction which is far too low especially at small x. A similar improvement was already noticed in the previous H1 three-jet analysis [16]. In the present analysis, extending to lower invariant three-jet masses, an excess is observed of the data compared to the O(α 3 s ) prediction at the lowest x ∼ 10 −4 . This excess is found to be enhanced and to become highly significant for topologies with two forward jets and one central jet. Excesses which are probably related, albeit less significant, were observed in the forward jet analyses from H1 [13] and ZEUS [15]. The new analysis corroborates the hypothesis that the DGLAP leading log Q 2 approximation starts to break down in the region of the lowest accessible x ∼ 10 −4 for Q 2 ≥ 4 GeV 2 , at least up to the order O(α 3 s ) for which calculations are presently available. In other words in a sizable fraction of events, which is much larger than predicted, two or more gluons are radiated from the initial state proton which are unordered in their transverse momentum, i.e. they all have relatively large transverse momenta. For events with at least four jets the O(α 3 s ) prediction is also too low, as expected, since the calculation can provide only a leading order estimate.
The new data presented here are also compared with two Monte Carlo simulation programs implementing hard QCD 2→2 processes complemented by parton showers modelling higher order effects. The RAPGAP [20] program, using k T ordered parton showers and including resolved photon processes, fails to describe the data. On the other hand the DJANGO [21] program with non k T -ordered gluon radiation as implemented in the colour dipole model (CDM) gives a remarkably good description of the measured cross sections for events with at least three and four jets and even for events with higher jet multiplicities. The remaining discrepancies at high p * T require further studies. The three-and four-jet production is investigated further by a detailed study of the jet topologies as represented by three-jet angles and scaled momenta. The best description for the case of three-jet production is obtained by the O(α 3 s ) NLOJET++ calculation, significantly better than that provided by CDM.       Figure 12: Differential cross sections for events with at least four jets as a function of the transverse momentum p * T 1 of the leading jet, the pseudorapidity difference of the leading and the fourth jet η 1 − η 4 (with p * T 1 > p * T 2 > p * T 3 > p * T 4 ), the scaled energy X ′ 2 and the angle θ ′ as defined in figure 3. For the determination of X ′ 2 and θ ′ the four jets are reduced to three by combining the two jets with the lowest dijet mass. The resulting three jets (labelled with ′′ ) are then sorted with respect to their energy in the three-jet centre of mass frame (E ′′ 1 > E ′′ 2 > E ′′ 3 ) and both variables are calculated as for the three-jet case, for instance . For the details of the plotting of the data see the caption of figure 9. The additional global normalisation error of the data ( +22 −19 %) is not displayed. The data are compared to the RAPGAP and CDM predictions. Both Monte Carlo predictions are normalised to the total cross section of the data. Tables   Cross sections for events with at least three  0.65 ± 0.04 Table 3: Bin-averaged differential cross sections for the selected events with at least three jets in the kinematic range listed in table 1. The cross sections are defined at the hadron level and are given as a function of the jet multiplicity N Jet and the Bjørken scaling variable x. The following cross section uncertainties are indicated: statistical (δ stat ), uncorrelated systematic δ uncorr and correlated systematic (δ corr ). The additional global normalisation uncertainty of +16 −14 % is not included in the table. The correction factors c had for the effect of hadronisation and the associated uncertainty are also given. They are applied to the NLOJET++ parton level calculations.

Cross sections for events with two forward jets and one central jet
x dσ/dx δ stat δ uncorr δ corr c had (pb) (%) (%) (%) [0.0001, 0.0002) 1.3 · 10 5 4 8 Cross sections for events with at least three jets and p * T 1 > 20 GeV  Table 9: Bin-averaged differential cross sections at the hadron level as a function of the Bjørken scaling variable x and the pseudorapidity η 1 of the leading jet in the lab frame, The cross sections are measured in the kinematic range listed in table 1. In addition the leading jet is required to have a large transverse momentum p * T 1 > 20 GeV in the γ * p centre of mass frame. The following cross section uncertainties are given: statistical (δ stat ), uncorrelated systematic δ uncorr and correlated systematic (δ corr ). The additional global normalisation uncertainty of +19 −14 % is not included in the table. The correction factors c had for the effect of hadronisation and the associated uncertainty are also given. They are applied to the NLOJET++ parton level calculations.