Measurement of charged particle multiplicity distributions in DIS at HERA and its implication to entanglement entropy of partons

Charged particle multiplicity distributions in positron-proton deep inelastic scattering at a centre-of-mass energy s=319\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{s}=319$$\end{document} GeV are measured. The data are collected with the H1 detector at HERA corresponding to an integrated luminosity of 136 pb-1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{-1}$$\end{document}. Charged particle multiplicities are measured as a function of photon virtuality Q2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Q^2$$\end{document}, inelasticity y and pseudorapidity η\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta $$\end{document} in the laboratory and the hadronic centre-of-mass frames. Predictions from different Monte Carlo models are compared to the data. The first and second moments of the multiplicity distributions are determined and the KNO scaling behaviour is investigated. The multiplicity distributions as a function of Q2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Q^2$$\end{document} and the Bjorken variable xbj\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_{\mathrm{bj}}$$\end{document} are converted to the hadron entropy Shadron\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$S_{\mathrm{hadron}}$$\end{document}, and predictions from a quantum entanglement model are tested.


Introduction
In the parton model [1][2][3] formulated by Bjorken, Feynman, and Gribov, the quarks and gluons of a nucleon are viewed as "quasi-free" particles probed by an external hard probe in the infinite momentum frame. The parton that participates in the hard interaction with the probe, e.g., a virtual photon, is expected to be causally disconnected from the rest of the nucleon. On the other hand, the parton and the rest of the nucleon have to form a colour-singlet state due to colour confinement. In order to further understand the role of colour confinement in high energy collisions, it has been suggested [4,5] that quantum entanglement of partons could be an important probe of the underlying mechanism. Measurements of charged particle multiplicities are proposed [5][6][7][8][9][10] to be related to the entanglement entropy predicted from the gluon density [5], as an indication of quantum entanglement of partons inside the proton.
Particle multiplicity distributions have previously been measured in deep inelastic scattering (DIS) at HERA [11][12][13][14][15]. In this paper, charged particle multiplicity distributions in positronproton (ep) DIS at √ s = 319 GeV are reported using high statistics data collected with the H1 detector. The phase space of the measurement of multiplicity distributions is defined in bins of the virtuality of the photon 5 < Q 2 < 100 GeV 2 , the inelasticity y variable 0.0375 < y < 0.6 and the pseudorapidity η of charged particles, |η lab | < 1.6 in the laboratory frame, and 0 < η * < 4 in the hadronic centre-of-mass (HCM) frame. 1 The first and the second moments of the multiplicity distributions, corresponding to the mean and the variance, are measured as a function of the hadronic centre-of-mass energy W , in different bins of Q 2 . The KNO scaling function [16] is also measured in different W and Q 2 regions.
The final-state hadron entropy, S hadron , as a function of the Bjorken variable x bj in different Q 2 bins is also measured. A relation between the S hadron and the initial-state parton entropy, S gluon , due to "parton liberation" [17] and "local parton-hadron duality (LPHD)" [18], is described by [5]: S hadron ≡ − ∑ P(N) ln P(N) = ln [xG(x, Q 2 )] ≡ S gluon .
(1) Here, P(N) is the charged particle multiplicity distribution measured in either the current fragmentation region or the target fragmentation region, where the gluon density xG(x, Q 2 ) is evaluated at x = x bj .

H1 detector
A full description of the H1 detector can be found elsewhere [22][23][24][25][26][27][28] and only the components most relevant for this analysis are briefly mentioned here. The coordinate system of H1 is defined such that the positive z axis is pointing in the proton beam direction (forward direction) and the nominal interaction point is located at z = 0. The polar angle θ is defined with respect to this axis. The pseudorapidity is defined to be η ≡ − ln (tan (θ /2)).
Charged particles are measured in the polar angle range 15°< θ < 165°using the central tracking detector (CTD), which is also used to reconstruct the interaction vertex. The CTD comprises two large cylindrical, concentric and coaxial jet chambers (CJCs), and the silicon vertex detector [25,26]. The CTD is operated inside a 1. 16 T solenoidal magnetic field. The CJCs are separated by a cylindrical drift chamber which improves the z coordinate reconstruction. A cylindrical multiwire proportional chamber [27], which is mainly used in the trigger, is situated inside the inner CJC. The trajectories of charged particles are measured with a transverse momentum resolution of σ (p T )/p T ≈ 0.2%/ GeV ⊕ 1.5%. The forward tracking detector (FTD) [28] measures the tracks of charged particles at polar angles 6°< θ < 25°. In the region of angular overlap, FTD and short CTD track segments are used to reconstruct combined tracks, extending the detector acceptance for well-reconstructed tracks. Both CTD tracks and combined tracks are linked to hits in the vertex detectors: the central silicon tracker (CST) [25,26], the backward silicon tracker (BST) and the forward silicon tracker (FST). These detectors provide precise spatial coordinate measurements and therefore significantly improve the primary vertex spatial resolution. The CST consists of two layers of double-sided silicon strip detectors surrounding the beam pipe covering an angular range of 30°< θ < 150°for tracks passing through both layers. The BST consists of six double wheels of silicon strip detectors measuring the transverse coordinates of charged particles. The FST design is similar to that of the BST and consists of five double wheels of single-sided silicon strip detectors. The lead-scintillating fibre calorimeter (SpaCal) [24] covering the region 153°< θ < 177.5°has electromagnetic and hadronic sections. The calorimeter is used to measure the scattered positron and the backward hadronic energy flow. The energy resolution for positrons in the electromagnetic section is σ (E)/E ≈ 7.1%/ E/ GeV ⊕ 1%, as determined in test beam measurements [29]. The SpaCal provides energy and time-of-flight information used for triggering purposes. A backward proportional chamber (BPC) in front of the SpaCal is used to improve the angular measurement of the scattered lepton. The liquid argon (LAr) calorimeter [30] covers the range 4°< θ < 154°a nd is used in this analysis in the reconstruction of the hadronic final state. It has an energy resolution of σ (E)/E ≈ 50%/ E/ GeV ⊕ 2% for hadronic showers, as obtained from test beam measurements [31].

Theoretical predictions
The DIS process is simulated by different Monte Carlo (MC) event generators, which include the hard scattering process and simulations of higher order QCD correction in form of parton showers and hadronisation. A brief description of the MC event generators is given below: • The RAPGAP 3.1 [19] MC event generator matches first order Quantum Chromodynamics (QCD) matrix elements to the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) [32][33][34][35] parton showers with strongly ordered transverse momenta of subsequently emitted partons. The factorisation and renormalisation scales are set to µ f = µ r = Q 2 +p 2 T , wherep T is the transverse momentum of the outgoing hard parton from the matrix element in the centre-of-mass frame of the hard subsystem. The CTEQ 6L [36] leading order parametrisation of the parton density function (PDF) is used.
• The DJANGOH 1.4 [20] MC event generator uses the Colour Dipole Model (CDM) as implemented in ARIADNE [37], which models first order QCD processes and creates dipoles between coloured partons. Gluon emission is treated as radiation from these dipoles, and new dipoles are formed from the emitted gluons from which further radiation is possible. The radiation pattern of the dipoles includes interference effects, thus modelling gluon coherence. The transverse momenta of the emitted partons are not ordered in transverse momentum with respect to rapidity, producing a configuration similar to the Balitsky-Fadin-Kuraev-Lipatov (BFKL) [38][39][40] treatment of parton evolution [41]. The CTEQ 6L [36] at leading order is used as the PDF.
• The PYTHIA 8 [21] MC event generator models the hard collision by LO-pQCD cross sections. In order to enable the PYTHIA 8 program to simulate DIS processes, the O(α 0 s ) process together with the DIrect REsummation (DIRE) parton shower [42] is applied. The parton shower recoil is treated such that the kinematic variables of the DIS process (y and Q 2 ) are unchanged. The PDF CTEQ 5L [36] is used.
DJANGOH and RAPGAP, as well as the PYTHIA 6 [43] photoproduction events, are also used together with the H1 detector simulation in order to determine acceptance, efficiency and backgrounds. as well as to estimate systematic uncertainties associated with the measurement. DJANGOH and RAPGAP are interfaced to HERACLES [44][45][46] to simulate the QED radiative effects. The generated events are passed through a detailed simulation of the H1 detector response based on the GEANT3 simulation program [47] and are processed using the same reconstruction and analysis chain as used for the data. For the determination of the detector effects both the RAPGAP and DJANGOH predictions are studied.

Event selection
The data set used for this analysis was collected with the H1 detector in the years 2006 and 2007 when positrons and protons were collided at energies of 27.6 GeV and 920 GeV, respectively. The integrated luminosity of the data set is 136 pb −1 [48]. DIS events were recorded using triggers based on electromagnetic energy deposits in the SpaCal calorimeter. The trigger inefficiency is determined using independently triggered data and is negligible in the kinematic region of the analysis.
The scattered positron is defined by the energy cluster in the SpaCal calorimeter with the highest transverse momentum. The energy of the cluster is further required to be larger than 12 GeV, and the radial position of the cluster is required to be between 15 and 70 cm with respect to the beam axis. The z coordinate of the event vertex is required to be within |z v | < 35 cm of the nominal interaction point. The E − p z variable is required to be between 35 and 75 GeV in order to reduce QED radiation and photoproduction backgrounds. Here E − p z is defined as the sum of E i − p z,i of both the scattered positron and the hadronic final-state (HFS) particles. The HFS particles are reconstructed using an energy flow algorithm [49][50][51]. This algorithm combines charged particle tracks and calorimetric energy clusters into hadronic objects, taking into account their respective resolution and geometric overlap, while avoiding double counting of energy.
Tracks measured in the CTD alone (central tracks) and the combination of CTD and FTD (combined tracks) are used in this analysis. Both central and combined tracks are required to have transverse momenta p T,lab > 150 MeV. Furthermore, the total momentum for the combined tracks is required to be larger than 500 MeV in order to ensure particles have enough momentum to cross the material between the CJC and FTD. The pseudorapidity of both central and combined tracks is required to be within |η lab | < 1.6. In addition, the central tracks are required to have their |dca sin θ | < 2 cm, where dca is the distance of closest approach with respect to the primary vertex and θ is the polar angle of the track. The innermost hit in the CTD is required for central tracks to be less than 50 cm away from the z axis, where the radial track length is required to be larger than 10 cm. Neutral particles are not considered in the multiplicity analysis. Using only tracks assigned to the primary event vertex, the contributions from in-flight decays of K 0 S , Λ, from photon conversions and from other secondary decays and interactions with detector material are reduced. Details of the track selection are described elsewhere [52].

Event reconstruction
The scattered positron and the HFS particles are used for the reconstruction of the kinematic variables, x bj , y and Q 2 . Similar to the H1 measurement of charged particle momentum spectra [52], the e-Σ method is used [53], where these variables are defined as Here, s is the ep centre-of-mass energy squared, E e is the incoming lepton energy, E e and θ e are the scattered positron energy and polar angle, respectively. The quantity Σ is defined as ∑ E i − p z,i , where the sum runs over all the HFS particles. This method provides an optimum in resolution of the kinematic variables and shows only little sensitivity to QED radiative effects. The hadronic centre-of-mass energy W is defined as: where M p is the proton rest mass.
The hadronic centre-of-mass frame is defined as the frame where p + q = 0, with p and q being the four-momenta of the proton and the virtual photon, respectively. In this frame, the positive z axis is aligned with the direction of the virtual photon.
In figure 1, distributions of the reconstructed quantities of Q 2 , y, η lab , and the number of charged particles N rec , are shown in comparison with predictions from DJANGOH and RAP-GAP. The panels of the Q 2 and y distribution also contain the contribution of expected photoproduction background estimated based on PYTHIA 6.4 [54]. The photoproduction background is found to be less than 0.5%, and therefore neglected in the subsequent analysis.

Experimental observables and kinematic phase space
For a given range in Q 2 and y, the probability P(N) is defined as the fraction of events for which N charged particles are produced in the specified η range relative to the total number of events in that Q 2 ,y range. Based on the distribution P(N), the first and second moments of the multiplicity distributions, the KNO function Ψ(z) and the final-state hadron entropy S hadron are defined as: Here, the sum runs over the number of charged particles and the z variable is equal to N/ N . These quantities are measured within the fiducial kinematic phase space listed in Table 1. The selection in the HCM frame relative to the laboratory frame differs only by the additional restriction of the charged particles properties to the current hemisphere, 0 < η * < 4.

Data correction
The detector-level charged track multiplicity distributions are corrected to stable particles with proper lifetime cτ > 10 mm including charged hyperons. The probabilities P(N) are derived from the distributions of events with N = 0, 1, 2, ... tracks reconstructed in the specified η range and with Q 2 and y reconstructed in the kinematic bin of interest. First, the observed event counts in the three-dimensional grid of N, Q 2 , y are unfolded, such that migrations between bins as well as efficiency and acceptance distortions are removed. The second step is to normalize the event counts as a function of N to the total number of events in each Q 2 , y bin and thus obtain the probabilities P(N). The last step is to correct for QED radiation from the electron line.
The unfolding is done within the TUnfold framework [55]. In order to better resolve migration effects, the number of bins in Q 2 , y and N is chosen to be higher when counting events in reconstructed quantities as compared to the truth quantities. The phase space in Q 2 and y is enlarged over the phase space given in Table 1, such that the measured phase space is guarded against migrations at the phase space boundaries. The unfolding turns out to be robust against variations of the regularisation scheme. The dominant systematic uncertainty in the detector correction procedure arises from constructing the matrix of migration with an alternative MC model, DJANGOH instead of the RAPGAP default.
In this paragraph QED radiative effects are discussed. In the radiative MC, the scattered positron at the particle level is defined to be a scattered positron with photons that are within a cone of 0.4 radian. The QED radiative effects are corrected for based on the MC event generator. The correction factors are derived bin-by-bin in each measured phase space at the particle level using radiative and non-radiative MC samples.
The binning of the multiplicity N at the particle level after the unfolding is made to be wider than at the detector level when N > 3, in order to avoid large negative bin-to-bin correlations. For those wide bins, the reported values of P(N) are defined as ∑ P(N)/∆, where ∆ is the number of distinct values of N included in the bin. In order to determine moments and the hadron entropy, the measured distribution is extrapolated within the wide bins using an exponentiated cubic spline. Correlations between the extrapolated values are taken into account. Bin centres of the wide bins are also reported.

Systematic uncertainties
Systematic uncertainties are studied based on systematic variations on fully corrected results, where details of the variations are listed below. The systematic uncertainties are found to depend on N, η lab , Q 2 and y. They are are observed to be similar in size for the HCM results as compared to those obtained in the laboratory frame. The following systematic uncertainty sources are considered in this analysis: • Sys.1 -Radiative corrections: the difference in correction factors between the MC generators, the DJANGOH and RAPGAP is taken as systematic uncertainty. It is found to be 1-2% for the P(N) distributions and up to 1% for their moments.
• Sys.2 -MC model dependence: the P(N) distributions and their moments are compared between results that are corrected by DJANGOH or RAPGAP. This leads to 1-4% systematic uncertainty on the multiplicity distributions and their moments.
• Sys.3 -SpaCal energy scale and angular resolution: a variation of 1% on the energy scale of the SpaCal [56] and 1 • on the angular direction is considered. The combined systematic uncertainty on the distribution of P(N) and its moments is found to be 1-5% and 2%, respectively.
• Sys.4 -Hadronic energy scale: a variation of 2% on the energy scale of hadronic finalstate objects [57] results in 1-3% systematic uncertainty on the multiplicity distribution and 2% on the moments.
• Sys.5 -Single track efficiency: 0.5% of tracking efficiency uncertainty on central tracks and 10% uncertainty on combined tracks are applied. This leads to 1-5% systematic uncertainty on the P(N) distributions, and 1% on the moments.
• Sys.6 -V0 particle contamination: 50% of uncertainty on tracks originating from K 0 s decay and 0.5% uncertainty on tracks from photon conversion and Dalitz-decays surviving the primary track selection result in 1-7% of systematic uncertainty on the P(N) distributions and up to 2% for the moments.
• Sys.7 -Diffractive contributions: variation of the diffractive MC contribution by a factor of 2 results in 1-5% systematic uncertainty on the P(N) distributions and up to 1% on the moments.
• Sys.8 -Extrapolation: values of P(N) that are not explicitly measured are extrapolated and 1% of difference is found for the mean, variance, and entropy with respect to the generator value based on MC.
Systematic uncertainties associated with photoproduction background are negligible. A summary of systematic uncertainties as a function of the multiplicity N can be found in Table 2.
The total systematic uncertainties are obtained by adding in quadrature all individual contributions. The tables of the results contain a complete breakdown of the systematic uncertainty contributions from different sources for each measured data point.  Table 2: Summary of the systematic uncertainties in this analysis.

Multiplicity distributions
The charged particle multiplicity distributions in ep DIS at √ s = 319 GeV are measured in different bins of Q 2 and y for particles with p T,lab > 150 MeV and |η lab | < 1.6 in the laboratory frame. The data are presented in figure 2, numerical values are reported in the appendix.
The data are compared with predictions from the DJANGOH, RAPGAP, and PYTHIA 8 event generators without simulation of QED radiation. The multiplicity distributions P(N) are found to broaden as y increases for fixed Q 2 . Since y can be related to the hadronic centre-of-mass energy W (cf. Eq. 3), the increase in multiplicity is qualitatively expected because more energy is available for hadronisation. On the other hand, the P(N) distributions are found to be almost independent of Q 2 for fixed y. Qualitatively, the MC predictions can describe the peak position of the multiplicity distributions well. However, they tend to underestimate the data both at low and high multiplicities, especially towards low Q 2 and high y. Among all MC models considered, the PYTHIA 8 model gives the poorest description and peaks significantly below the data, especially at high y.
In Figures. 3 to 6 the charged particle multiplicity distributions P(N) are presented in four bins of Q 2 . In each Figure, the P(N) distributions are shown differentially in bins of y (identical binning as in figure 2) and in three different ranges of η lab . Numerical values are given in the appendix. The overlapping η lab ranges are chosen such that a pseudorapidity window of 1.4 units around the scattered parton direction in the leading order Quark Parton Model (QPM) picture can be selected. Similar to the measurements over the full η lab range, the P(N) distributions are found to broaden as y increases, independent of the considered η lab ranges. For fixed y, the P(N) distributions also broaden as η lab increases. Similar to the results presented in figure 2, the MC models underestimate the high multiplicity tail, where the deviation is found to be the strongest at low Q 2 .
In figure 7, the multiplicity distributions P(N) are presented with the additional restriction in the HCM frame 0 < η * < 4, mainly selecting particles originating from the current hemisphere. Numerical values are given in the appendix. Predictions obtained from DJANGOH, RAPGAP, and PYTHIA 8 are compared with data. Qualitatively, the comparisons are very similar to those of Figure. 2. For the results presented in the subsequent sections, only RAPGAP is compared to the data.

Moments of multiplicity distributions
The mean multiplicity N as a function of W is presented in figure 8 and tables 3 to 4, for both the full phase space and with an additional restriction in η * (see Table 1). Predictions from the RAPGAP model are also shown. The average multiplicity N rises with the hadronic centreof-mass energy W , which is in agreement with previous observations [12,13,15,[58][59][60][61]. For higher Q 2 , the increase in W is faster than that at lower Q 2 . For high W , the Q 2 dependence is observed to be stronger in the full phase space than with the η * restriction of the current hemisphere. The RAPGAP prediction yields a reasonable description of the data. Only at high W and low Q 2 the MC tends to underestimate the data.
Similarly, in figure 9 and tables 5 to 6, second moments of the multiplicity distributions (variances) are presented as a function of W , for both the full phase space and with the additional restriction in η * . The variances rise strongly with W , almost independent of Q 2 within uncertainties. The restriction of the current hemisphere has little inpact on the variance. For Q 2 > 40 GeV 2 , RAPGAP describes the data reasonably well. However, towords high W , the MC not only underestimates the data but also shows a Q 2 dependence, which is absent in data. This effect is more pronounced when restricting the analysis to the current hemisphere.

The KNO-scaling function Ψ(z)
In order to further study the multiplicity distribution, the KNO function Ψ(z) measured as a function of z = N/ N is shown in different bins of Q 2 in figure 10. Numerical values are given in the appendix. The analysis is restricted to the current hemisphere (0 < η * < 4). In all Q 2 bins, KNO scaling is observed, in broad agreement with many past experiments. However, in proton-proton collisions at LHC energies, violations of KNO scaling have been reported recently [60].

Entropy
It was recently suggested [5,10] that the final-state hadron entropy S hadron calculated from charged particle multiplicity distributions might be related to the entanglement entropy of gluons S gluon at low x bj (Eq. 1). In figure 11 and Table 7, S hadron is studied as a function of x bj in different Q 2 bins. Similar to the observable studied in Ref. [10], a moving η lab window depending on x bj is chosen to match the rapidity of the scattered quark in a leading order QPM picture. The respective selected pseudorapidity window in the laboratory frame is indicated in Table 7. The hadron entropy observed in data is cosistent with being constant in x bj , increasing only slightly with Q 2 . The same x bj and Q 2 dependence is also observed for RAPGAP, however with slightly smaller values of S hadron . Predictions for the entanglement entropy S gluon based on the gluon density xG(x, Q 2 ) are also shown at various Q 2 values, which correspond to the lower boundaries of the Q 2 bins in data. Here the PDF set HERAPDF 2.0 [62] at leading order is used. Neither the dependence on x bj nor the magnitude agrees with data. Thus the prediction S hadron = S gluon (Eq. 1) is not confirmed by the present measurement.
To investigate further, the hadron entropy determined in the current hemisphere, 0 < η * < 4, is presented in figure 12 and Table 8. The hadron entropy based on multiplicity distributions is shown as a function of x bj in different bins of Q 2 . In contrast to Figure 11, there is no moving pseudorapidity range applied. Independent of Q 2 , the measured hadron entropy, S hadron , rises with decreasing x bj with similar slopes in different Q 2 bins. Also shown are the predictions from RAPGAP which closely follow the data at high Q 2 but underestimate the data at low Q 2 . Since S hadron is calculated from the multiplicity distributions, this behaviour of the MC is related to what has been observed for the moments, figure 8 and 9. The predictions from the entanglement approach based on the gluon density again fail to describe S hadron in magnitude. However, at low Q 2 the slope of S gluon has some similarities with that observed for S hadron , while it becomes steeper than observed with increasing Q 2 .

Summary
The charged particle multiplicity distributions P(N) are measured in deep inelastic scattering at √ s = 319 GeV using the H1 detector at HERA. The integrated luminosity used in this analysis is 136 pb −1 , recorded in the years 2006 and 2007 in positrons scattering off protons. The P(N) distributions are measured in bins of Q 2 , y and pseudorapidity η. Predictions based on simulations of partonic tree level matrix elements, supplemented by parton shower and hadronisation are generally found to be consistent with the measurement at low multiplicities while they underestimate the data in the high multiplicity regions, especially at low Q 2 . For measurements of moments, the predictions generally describe the data well at low hadronic centre-of-mass energy W and high Q 2 but less so at high W and low Q 2 . In addition, KNO scaling is observed within the kinematic phase space of the analysis.
The measurement of the charged particle multiplicity distributions is also used to test for the first time predictions based on quantum entanglement on sub-nucleonic scales in deep-inelastic ep scattering. The predictions from the entropy of gluons are found to grossly disagree with the hadron entropy obtained from the multiplicity measurements presented here, and therefore the data do not support the basic concept of equality of the parton and hadron entropy with the current level of theory development.
The measurements reported in this paper not only provide valuable information for better understanding particle production mechanisms, but also set an important testing ground for the development of new concepts, like quantum entanglement at sub-nucleonic scales.  Figure 1: Reconstructed momentum transfer Q 2 , inelasticity y, pseudorapidity η lab , and charged particle multiplicity N rec for data (open circle), and the DJANGOH (red line) and the RAPGAP (blue line) MC models. The phase space restrictions are given in Table 1. The photoproduction background simulated from PYTHIA 6.4 [54] is shown for the Q 2 and y distributions. Error bars indicate the statistical uncertainty of data.  Table 1. Predictions from DJANGOH, RAPGAP and PYTHIA 8 are also shown. The total uncertainty is denoted by the error bars.  Table 1. Different panels correspond to different η lab and y bins, as indicated by the text in the figure. Predictions from DJANGOH, RAPGAP and PYTHIA 8 are also shown. The total uncertainty is denoted by the error bars.  Table 1. Different panels correspond to different η lab and y bins, as indicated by the text in the figure. Predictions from DJANGOH, RAPGAP and PYTHIA 8 are also shown. The total uncertainty is denoted by the error bars.    Table 1. Different panels correspond to different Q 2 and y bins, as indicated by the text in the figure. Predictions from DJANGOH, RAPGAP and PYTHIA 8 are also shown. The total uncertainty is denoted by the error bars.  Table 1. The corresponding y is indicated by the top axis for each measured W . Predictions from the RAPGAP model are shown by dashed lines. The total uncertainty is represented by the error bar.     Table 1. The total uncertainty is denoted by the error bars. For each x bj , the multiplicity is determined in a dedicated pseudorapidity window as discussed in the text. Further phase space restrictions are given in Table 1. Predictions for S hadron from the RAPGAP model and for the entanglement entropy S gluon based on an entanglement model are shown by the dashed lines and solid lines, respectively. For each Q 2 range, the value of the lower boundary is used for predicting S gluon . The total uncertainty on the data is represented by the error bars. Here, a restriction to the current hemisphere 0 < η * < 4 is applied. Further phase space restrictions are given in Table 1. Predictions for S hadron from the RAPGAP model and for the entanglement entropy S gluon based on an entanglement model are shown by the dashed lines and solid lines, respectively. For each Q 2 range, the value of the lower boundary is used for predicting S gluon . The total uncertainty on the data is represented by the error bars.       Table 8: Hadron entropy derived from multiplicity distributions measured in ep DIS at √ s = 319 GeV as a function of x bj . The measurement is repeated in four ranges of Q 2 as indicated. The phase-space is further restricted as shown in table 1. For this dataset, the entropy is determined with the track pseudorapidities in the hadronic centre-of-mass frame are restricted to the range 0 < η * < 4

A Charged particle multiplicities and KNO function
Numerical values of the charged particle multiplicity measured in the pseudorapidity range −1.6 < η lab < 1.6 are presented in four ranges of Q 2 and are further subdivided into four y ranges. The four Q 2 ranges are shown in the tables 9, 10, 11, and 12, respectively. Further phase space restrictions are discussed in table 1. These restrictions apply to all measurements presented in the following.
The measurement is repeated in three overlapping subranges of η lab . There are a total of twelve tables, each corresponding to selected range in both Q 2 and η lab . The range 5 < Q 2 < 10 GeV 2 is given in tables 13, 14, and 15. The range 10 < Q 2 < 20 GeV 2 is given in tables 16, 17, and 18. The range 20 < Q 2 < 40 GeV 2 is given in tables 19, 20, and 21. The range 40 < Q 2 < 100 GeV 2 is given in tables 22, 23, and 24.
The measurement is also repeated over the full available pseudorapidity range −1.6 < η lab < 1.6 in the laboratory frame, however with the extra restriction of the track selection to the current hemisphere 0 < η * < 4. Again the data are given in four Q 2 ranges, shown in the tables 25, 26, 27, and 28, respectively. For these data, the KNO function Ψ(z) is also determined. It is reported in the tables 29, 30, 31, and 32.
sys  Table 9: Charged particle multiplicity distributions P(N) measured in ep DIS at √ s = 319 GeV with in the ranges of photon virtuality 5 < Q 2 < 10 GeV 2 and charged particle pseudorapidity range |η lab | < 1.6, presented in four ranges of inelasticity y as indicated. The phasespace is further restricted as shown in  Table 10: Charged particle multiplicity distributions P(N) measured in ep DIS at √ s = 319 GeV with in the ranges of photon virtuality 10 < Q 2 < 20 GeV 2 and charged particle pseudorapidity range |η lab | < 1.6, presented in four ranges of inelasticity y as indicated. The phase-space is further restricted as shown in table 1.  Table 11: Charged particle multiplicity distributions P(N) measured in ep DIS at √ s = 319 GeV with in the ranges of photon virtuality 20 < Q 2 < 40 GeV 2 and charged particle pseudorapidity range |η lab | < 1.6, presented in four ranges of inelasticity y as indicated. The phase-space is further restricted as shown in table 1.  Table 12: Charged particle multiplicity distributions P(N) measured in ep DIS at √ s = 319 GeV with in the ranges of photon virtuality 40 < Q 2 < 100 GeV 2 and charged particle pseudorapidity range |η lab | < 1.6, presented in four ranges of inelasticity y as indicated. The phase-space is further restricted as shown in table 1.  Table 13: Charged particle multiplicity distributions P(N) measured in ep DIS at √ s = 319 GeV with in the ranges of photon virtuality 5 < Q 2 < 10 GeV 2 and charged particle pseudorapidity range −1.2 < η lab < 0.2, presented in four ranges of inelasticity y as indicated. The phase-space is further restricted as shown in table 1.
sys  Table 14: Charged particle multiplicity distributions P(N) measured in ep DIS at √ s = 319 GeV with in the ranges of photon virtuality 5 < Q 2 < 10 GeV 2 and charged particle pseudorapidity range −0.5 < η lab < 0.9, presented in four ranges of inelasticity y as indicated. The phase-space is further restricted as shown in table 1.
sys  Table 15: Charged particle multiplicity distributions P(N) measured in ep DIS at √ s = 319 GeV with in the ranges of photon virtuality 5 < Q 2 < 10 GeV 2 and charged particle pseudorapidity range 0.2 < η lab < 1.6, presented in four ranges of inelasticity y as indicated. The phase-space is further restricted as shown in  Table 16: Charged particle multiplicity distributions P(N) measured in ep DIS at √ s = 319 GeV with in the ranges of photon virtuality 10 < Q 2 < 20 GeV 2 and charged particle pseudorapidity range −1.2 < η lab < 0.2, presented in four ranges of inelasticity y as indicated. The phase-space is further restricted as shown in  Table 17: Charged particle multiplicity distributions P(N) measured in ep DIS at √ s = 319 GeV with in the ranges of photon virtuality 10 < Q 2 < 20 GeV 2 and charged particle pseudorapidity range −0.5 < η lab < 0.9, presented in four ranges of inelasticity y as indicated. The phase-space is further restricted as shown in  Table 18: Charged particle multiplicity distributions P(N) measured in ep DIS at √ s = 319 GeV with in the ranges of photon virtuality 10 < Q 2 < 20 GeV 2 and charged particle pseudorapidity range 0.2 < η lab < 1.6, presented in four ranges of inelasticity y as indicated. The phase-space is further restricted as shown in  Table 19: Charged particle multiplicity distributions P(N) measured in ep DIS at √ s = 319 GeV with in the ranges of photon virtuality 20 < Q 2 < 40 GeV 2 and charged particle pseudorapidity range −1.2 < η lab < 0.2, presented in four ranges of inelasticity y as indicated. The phase-space is further restricted as shown in  Table 20: Charged particle multiplicity distributions P(N) measured in ep DIS at √ s = 319 GeV with in the ranges of photon virtuality 20 < Q 2 < 40 GeV 2 and charged particle pseudorapidity range −0.5 < η lab < 0.9, presented in four ranges of inelasticity y as indicated. The phase-space is further restricted as shown in table 1.  Table 21: Charged particle multiplicity distributions P(N) measured in ep DIS at √ s = 319 GeV with in the ranges of photon virtuality 20 < Q 2 < 40 GeV 2 and charged particle pseudorapidity range 0.2 < η lab < 1.6, presented in four ranges of inelasticity y as indicated. The phase-space is further restricted as shown in  Table 22: Charged particle multiplicity distributions P(N) measured in ep DIS at √ s = 319 GeV with in the ranges of photon virtuality 40 < Q 2 < 100 GeV 2 and charged particle pseudorapidity range −1.2 < η lab < 0.2, presented in four ranges of inelasticity y as indicated. The phase-space is further restricted as shown in  Table 23: Charged particle multiplicity distributions P(N) measured in ep DIS at √ s = 319 GeV with in the ranges of photon virtuality 40 < Q 2 < 100 GeV 2 and charged particle pseudorapidity range −0.5 < η lab < 0.9, presented in four ranges of inelasticity y as indicated. The phase-space is further restricted as shown in  Table 24: Charged particle multiplicity distributions P(N) measured in ep DIS at √ s = 319 GeV with in the ranges of photon virtuality 40 < Q 2 < 100 GeV 2 and charged particle pseudorapidity range 0.2 < η lab < 1.6, presented in four ranges of inelasticity y as indicated. The phase-space is further restricted as shown in table 1. 5 < Q 2 < 10 GeV 2 and 0 < η * < 4 0.0375 < y < 0.075 stat.

Multiplicity distributions in DIS
sys  Table 25: Charged particle multiplicity distributions P(N) measured in ep DIS at √ s = 319 GeV with in the ranges of photon virtuality 5 < Q 2 < 10 GeV 2 and charged particle pseudorapidity range 0 < η * < 4, presented in four ranges of inelasticity y as indicated. The phasespace is further restricted as shown in  Table 26: Charged particle multiplicity distributions P(N) measured in ep DIS at √ s = 319 GeV with in the ranges of photon virtuality 10 < Q 2 < 20 GeV 2 and charged particle pseudorapidity range 0 < η * < 4, presented in four ranges of inelasticity y as indicated. The phase-space is further restricted as shown in table 1.
sys  Table 27: Charged particle multiplicity distributions P(N) measured in ep DIS at √ s = 319 GeV with in the ranges of photon virtuality 20 < Q 2 < 40 GeV 2 and charged particle pseudorapidity range 0 < η * < 4, presented in four ranges of inelasticity y as indicated. The phase-space is further restricted as shown in table 1.
sys  Table 28: Charged particle multiplicity distributions P(N) measured in ep DIS at √ s = 319 GeV with in the ranges of photon virtuality 40 < Q 2 < 100 GeV 2 and charged particle pseudorapidity range 0 < η * < 4, presented in four ranges of inelasticity y as indicated. The phase-space is further restricted as shown in table 1. 3 Table 29: KNO function Ψ(z) of charged particle multiplicity distributions measured in ep DIS at √ s = 319 GeV with in the ranges of photon virtuality 5 < Q 2 < 10 GeV 2 and charged particle pseudorapidity range 0 < η * < 4, presented in four ranges of inelasticity y as indicated. The phase-space is further restricted as shown in table 1.