High energy QCD: multiplicity distribution and entanglement entropy

In this paper we show that QCD at high energies leads to the multiplicity distribution $\frac{\sigma_n}{\sigma_{ \rm in}}\,\,=\,\,\frac{1}{N}\,\Lb \frac{N\,-\,1}{N}\Rb^{n - 1}$, (where $N$ denotes the average number of particles), and to entanglement entropy $S \,=\,\ln N$, confirming that the partonic stat at high energy is maximally entangled. However, the value of $N$ depends on the kinematics of the parton cascade. In particular, for DIS$N = xG(x,Q)$ , where $xG$ is the gluon structure function, whil for hadron-hadron collisions, $N \propto Q^2_S(Y)$, where $Q_s$ denotes the saturation scale. We checked that this multiplicity distribution describes the LHC data for low multiplicities $n \,<\,(3 \div 5)\,N$, exceeding it for larger values of $n$. We view this as a result of our assumption, that the system of partons in hadron-hadron collisions atc.m. rapidity $Y=0$ is dilute. We show that the data can be described at large multiplicities in the parton model, if we do not make this assumption.


arXiv:2006.11793v1 [hep-ph] 21 Jun 2020
This relation materialized as the resolution of the following difficulty in our understanding of high energy scattering: on one hand, the proton is a pure state and it is described by a completely coherent wave function with zero entropy, but, on the other hand, the DIS experiments are successfully described, treating the proton as a incoherent collection of quasi-free partons. This ensemble has non vanishing entropy, and Ref. [5] proposes, that the origin of this entropy is the entanglement between the degrees of freedom one observes in DIS (partons in the small spatial region of the proton), and the rest of the proton wave function, which is not measured in the DIS experiments.
In other words, the hadron in the rest frame is described by a pure state |ψ with density matrixρ = |ψ ψ| and zero von Neumann entropy S = −tr [ρ lnρ] = 0. In DIS at Bjorken x and momentum transfer q 2 = −Q 2 probes only a part of the proton's wave function; let us denote it A. In the proton's rest frame the DIS probes the spatial region A localized within a tube of radius ∼ 1/Q and length ∼ 1/(mx) [22,23], where m is the proton's mass. The inclusive DIS measurement thus sums over the unobserved part of the wave function localized in the region B complementary to A, so we have access only to the reduced density matrixρ A = tr Bρ , and not the entire density matrixρ = |ψ ψ|. In Ref. [5] it is proposed that Eq. (1), in spite of its general form, means that we can estimate the entropy and multiplicity distribution of the produced gluons using the parton wave function in the initial state. In addition , we can obtain a thermal distribution of the produced particles in the high energy collision in spite of the fact, that the number of secondary interactions in proton-proton collisions is rather low, and cannot provide the thermalization due to the interaction in the final state. It has been demonstrated in Refs. [6,10,11,16] that these ideas are in qualitative and, partly, in quantitative agreement with the available experimental data.
We believe that the CGC/saturation approach provides a more general pattern [43,44] for the treatment of high energy QCD. However, in this paper we restrict ourself to the BFKL Pomeron calculus, which has a more direct correspondence with the parton approach, and has been used in Ref. [5].
Fortunately, in Ref. [44] it was shown, that these two approaches are equivalent for the description of the scattering amplitude where ∆ BFKL denotes the intercept of the BFKL Pomeron.
The main difference between the CGC approach and the parton QCD cascade for the topics dealt with in this paper, is the fact that the CGC approach generates the non-diagonal elements of the density matrix (see for example Refs. [3,9,12,17]), while for the parton cascade and, generally, in the BFKL Pomeron calculus the density matrix is diagonal. Since in DIS experiments we can only measure the diagonal elements of density matrix, we can introduce in the framework another kind of entropy:" the entropy of ignorance" [17], which characterized this lack of knowledge of the actual density matrix in DIS experiments. We will show below that the McLerran-Venugopalan approach [25] ,that is used in Ref. [17], leads to the same multiplicity distribution as the parton cascade.
The paper is organized as follows: In the next section we consider the entropy and multiplicity distributions in the QCD parton cascade. We show that in spite of the fact that in different kinematic regions, the QCD cascade leads to a different energy and dipole size dependence of the mean multiplicity, and the multiplicity distribution has a general form: where N is the average number of partons. The entanglement entropy is equal to S parton cascade = ln N , confirming that the partonic state at high energy is maximally entangled [5]. In the case of DIS we argue that N is equal to the gluon structure function, but we can only prove the multiplicity distribution of Eq. (3), for the BFKL evolution of this structure function. In section III we show that the CGC approach leads to the multiplicity distribution of Eq. (3). In section IV we consider hadron-hadron scattering. In the range of energy given by Eq. (2) we use the Mueller-Patel-Salam-Iancu(MPSI) [36,49] approach, using the formalism of Ref. [41]. We show that in the framework of this approach we have the distribution of Eq. (3) which can describe the experimental data for sufficiently low multiplicities n ≤ (3 ÷ 5) < n > . However, we fail to describe the data for larger n. We conclude that the main assumption of the MPSI approach, that a system of dilute partons are produced in the c.m. rapidity Y=0, is not valid for large multiplicities at high energies of the LHC. Unfortunately, at the moment we have no theoretical tool to treat this scattering. However, in Ref. [50] an approach has been suggested, which allows us to describe the dense system of partons in hadron-hadron collisions, as well as the dilute one. Developing this approach for the multiplicity distribution, we are able to describe the data for large multiplicities. We summarize our results in the conclusions.

II. THE QCD PARTON CASCADE
A. QCD cascade for fast moving large dipole

General approach
As discussed in Refs. [24,26,40,41] the parton cascade can be written in the following form (see Fig. 1): is the probability to have n-dipoles of size r i , at impact parameter b i and at rapidity is a typical cascade equation in which the first term describes the reduction of the probability to find n dipoles due to the possibility that one of n dipoles can decay into two dipoles of arbitrary sizes , while the second term, describes the growth due to the splitting of (n − 1) dipoles into n dipoles. The initial condition for the DIS scattering is which corresponds to the fact that we are discussing a dipole of definite size which develops the parton cascade. Since P n (Y ; {r i }) is the probability to find dipoles {r i }, we have the following sum rule i.e. the sum of all probabilities is equal to 1. This QCD cascade leads to Balitsky-Kovchegov (BK) equation [24,27,28] for the amplitude and gives the theoretical description of the DIS. We introduce the generating functional [26] 1 In the lab. frame rapidity Y is equal to Y = y dipole r − y dipoles r i , where y dipole r is the rapidity of the incoming fast dipole and y dipole r i is the rapidity of dipoles r i .
where u (r i b i ) ≡ = u i is an arbitrary function. The initial conditions of Eq. (5) and the sum rules of Eq. (6) take to following form for the functional Z: Multiplying both parts of Eq. (4) by n i=1 u (r i b i ) and integrating over r i and b i we obtain the following linear functional equation [41]; Searching for the solution of the form Z ([u(r i , b i , Y )]) for the initial conditions of Eq. (8a), Eq. (9a) can be re-written as the non-linear equation [26]: (10) Therefore, the QCD parton cascade of Eq. (4) takes into account non-linear evolution. However, to obtain the BK equation for the scattering amplitude we need to introduce the scattering amplitude γ (r i , b), for the interaction of the dipole with the target at low energies. Using these amplitudes we can obtain the non-linear BK equation from Eq. (10), since [28] Using Eq. (9a) and Eq. (11) we derive the BK equation in the standard form:

Several first iterations
Our goal is to find the solution to Eq. (4). In particular, for the multiplicity distribution and for the entropy we wish to find P n is the probability to find n dipoles of all possible sizes at the same values of the impact parameters and, being such, it gives σ n /σ in , which is the multiplicity distribution in the QCD parton cascade. The initial and boundary conditions for P n (Y, r) follows from Eq. (5) and Eq. (6) and take the form: First, let us find P 1 (Y, r). The equation for P 1 has the form: with the initial condition Therefore for P 1 (Y, r, b) the equation takes the following form: with the solution: The equation for P 2 (Y, r, b, r 1 , b , r 2 , b ) has the following form: First, let us estimate the value of ω G (r) which is given by Eq. (9b): Hence, only dipoles of size smaller than r, contribute to the value of ω G (r).
We suggest that the solution to Eq. (19) has the following form: where Θ (z) denotes the step function:Θ (z) = 1 for z > 0, and Θ (z) = 0 for z < 0. For the solution of Eq. (21) we can obtain the equation for P 2 (Y, r, b), integrating both parts of Eq. (19) over b ,r 1 and r 2 . It has the form: Using Eq. (18) we obtain that One can see that Eq. (23) gives P 2 (Y = 0, r) = 0 in accord with Eq. (14). For small ω G (r) Y 1 in the parton cascade only two terms exist: P 1 (Y, r) and P 2 (Y, r), and Eq. (6) reduces to Eq. (24) shows that P n (Y, r, b, {r i , b }) are negligibly small for dipoles with large sizes ( r i > r).

Solution
We suggest to look for the general solution in the form: For such a solution we can obtain from Eq. (4) the following equations for P n (Y, r): Introducing the Laplace transform: we re-write Eq. (26) in the form: Eq. (28) has the solution: Taking the inverse Laplace transform of function Hence we have the following solution: It is easy to see that Eq. (31) satisfies the initial conditions and the sum rules of Eq. (14).

B. Multiplicity distribution and entropy of the parton cascade
As has been mentioned and, therefore, determines the multiplicity distribution. Calculating the average N we see that this distribution can be written in the form: where we have denotedN = N − 1.
We can compare Eq. (34) with the general form of the negative binomial distribution (NBD) One can see that Eq. (34) can be re-written as where σ n is the cross section for producing n hadrons in a collision, and σ in is the inelastic cross section. Therefore at largeN our distribution is close to the negative binomial distribution with number of failures r = 1 and with probability of success p =N / N + 1 . Eq. (34) coincides with the multiplicity distribution in the parton model (see below and Ref. [5]). The difference is only in the expression for the average multiplicity (N ). Having Eq. (34) we can calculate the von Neumann entropy of the parton cascade (see Eq. (1)), given by the Gibbs formula where p n is the probabilities of micro-states. In the parton cascade we can identify p n with P n (Y, r) , reducing Eq. (38) to the following expression: Eq. (38) shows that at large Y all probabilities P n become equal and small, of the order of P n ∼ 1 N . It is well known that this equipartitioning of micro-states maximizes the von Neumann entropy and describes the maximally entangled state. We thus conclude that at large Y the fast dipole represents a maximally entangled quantum state of partons.

C. QCD motivated parton model
The main assumption of the parton model [19][20][21] is that all partons have average transverse momentum which does not depend on energy. Therefore, we can obtain the parton model from the QCD cascade assuming that the unknown confinement of gluons leads to the QCD cascade for the dipole of fixed size. In this case the cascade equation take the following form: where P n (Y ) is the probability to find n dipoles (of a fixed size in our model) at rapidity Y and ∆ = ω G (r = r 0 ).
Using the Laplace transform of Eq. (27) we obtain the solution to Eq. (39) in the form: Using Eq. (30), we see that the solution of Eq. (39) has the form: which is a direct generalization of Eq. (31). In Eq. (41), N is the average number of the partons, which is equal to N = exp (∆ Y ). In the parton model the average number of partons is related to the deep inelastic structure function. Therefore, in Ref. [5] it is assumed that where xG is the gluon structure function. In this case ∆ can be identified with the intercept of the BFKL Pomeron [31]. One can see that N = exp (ω G (r) Y ) this is certainly not the same as a solution of the QCD evolution for the gluon structure function.

D. Mutiplicity distribution for the parton cascade in DIS
As we have seen (see Eq. (25)) our solution describes the evolution in the system of partons with smaller sizes of dipoles than the initial fastest dipole. On the other side, only partons of larger than initial parton size, contribute to the structure function. In double log approximation the emission of such dipoles leads to xG ∝ exp 2 ᾱ S Y ln (R 2 /r 2 ) , where R is the size of the target. We need to consider how the produced dipoles interact with the target. We measure the gluon structure function in the experiment in which r 2 of the fastest dipole is about r 2 sim1/Q 2 R 2 . Since all r i ≤ r the n produced dipoles interact with the target at rapidity Y and since r i < R the amplitude of this interaction is proportional to r 2 i /R 2 . This fact completely changes the structure of the cascade. Let us illustrate this, considering P 2 (Y, r; r 1 , r 2 ). The amplitude of interaction is proportional to The term with gluon ω G (r i ) does not contribute since as we have discussed only r i > r contribute in this term. Therefore the equation reduces to the following one: Therefore, we infer that dipole sizes larger than r, contribute to the scattering amplitude of our interest, and lead to large ρ 2 .
Generally speaking the scattering amplitude can be written in the form [28,41]: where N (Y 0 , r i , b i ) is the amplitude of the interaction of dipole r i with the target at low energy Y = Y 0 , and the n-dipole densities in the projectile ρ p n (r 1 , b 1 , . . . , r n , b n ) are defined as follows: For ρ n we obtain [41] : For ρ 1 we have the linear equation: The physical meaning of ρ p 1 is clear from Eq. (45): it is the mean number of dipoles with size r 1 that have been produced. The multiplicity, which we needed in DIS, is the number of dipoles with sizes larger than r ∼ 1/Q. It is equal to where ξ = ln(1/r 2 ).
In the double log approximation (DLA) of perturbative QCD, Eq. (48) can be re-written in the form It is worthwhile mentioning that N p 1 is the gluon structure function in the DLA. Equation forρ p 2 (Y ; r 1 , r 2 , b) has the form 2 : However, to find the multiplicity distribution we need to introduce moments (see Eq. (13)): N p 2 is equal to which gives n(n−1)
Equation for ρ p 2 can be re-written in the following form forρ p 2 : or in DLA it takes the form: Note, that the gluon reggeization does not contribute in DLA, since it describes the contribution of distances r i < r (see discussion above). Integrating Eq. (55) over ξ 1 and ξ 2 we obtain: The general solution to Eq. (56) has a form: is the solution of the homogenous equation: The solution of Eq. (57) has the form: where ω(γ) =ᾱ S γ is the DLA limit of the BFKL kernel: where ψ(z) is Euler gamma function (see [48] formula 8.36). We select n in (γ) = 1/γ, since at Y = 0 we have only one dipole and N 2 = < n(n − 1)/2 >= 0. Taking the integral over γ using the method of steepest descent, we obtain First we wish to note that Eq. (60) leads to N p,homog However, in the diffusion approximation for the BFKL kernel for N P n . Therefore, for large Y as well as for large n we, indeed, can consider ν SP → 0. For this special case we have Comparing with the multiplicity distribution of Eq. (34), one can see that Eq. (57) gives the factorial moment n(n−1) 2 of this distribution with < n >= N p 1 . For N p n , the equation follows from Eq. (53) which in DLA takes the form: The solution to this equation has the form: with ω (γ) =ᾱ S /γ in DLA. Comparing Eq. (65) with the moments M q = < n! q!(n−q)! > for the multiplicity distribution of Eq. (34): one can see that N p n (Y, ξ) coincide with these moments only if we take into account the main exponential behaviour at γ = 1/2. As we have seen above, for large n the value of the saddle point for ν (see Eq. (62) ) indeed approaching zero.
Hence, we infer that the QCD parton cascade in DIS leads to the multiplicity distribution of Eq. (34) with N = xG (Q, x) at x → 0, as is expected from the parton model of section II-E, but xG should satisfy the BFKL evolution equation, and the accuracy of Eq. (34) is not very precise at small n.

III. MULTIPLICITY DISTRIBUTION IN COLOUR GLASS CONDENSATE (CGC) APPROACH
In Ref. [17] the density matrix is calculated in the CGC approach, using the CGC wave function from Refs. [46,47]. In CGC approach the large fraction of momentum is carried by the valence quarks and gluons. These fast partons emit low energy gluons whose lifetime is much shorter than the valence partons. In other words the valence ("hard") partons can be treated as static sources of the soft gluons. The wave function of such a system of partons can be written in the form: where |v > characterizes the valence degrees of freedom, while |s > denotes the wave function of the soft gluon in the presence of the valence partons. Sign ⊗ does not denote a direct product, since the wave function of a soft gluon depends on the valence degrees of freedom. Using that where φ i (k) ≡ a + i + a(−k) and b i a = gρ a (k) iki k 2 + ..( ρ a is the charge density of the valence partons), and McLerran-Venugopalan (MV) model for wave function |v >, the density matrix: is calculated in Ref. [17]. The result of these calculation is the following: where g is QCD coupling constant, and µ 2 determines the colour charge density in the valence wave function in the MV model [25].
For the multiplicity distribution we only need the diagonal elements of the density matrix with l = α and m = β, and the multiplicity is n = l + m. Plugging these l and m into Eq. (70) we obtain for the multiplicity distribution Calculating average n = N we obtain and the multiplicity distribution can be re-written in the form: We stress that Eq. (74) coincides with Eq. (34), that we derived for the QCD parton cascade.

A. The interaction of two dipoles at high energies
We first consider the high energy interactions of two dipoles with sizes r and R and with r ∼ R. In Ref. [44] it is shown that in the limited range of rapidities, which is given by Eq. (1), we can safely apply the Muller, Patel, Salam and Iancu approach for this scattering [36,49](see Fig. 3-a). The scattering amplitude in this approach can be written FIG. 3: Scattering amplitude for the interaction of two dipoles with sizes: r and R at high energy in MPSI approach (see Fig. 3a and Fig. 3-b). The amplitudes of interaction of two dipoles in the Born approximation of perturbative QCD ( N (ri, r i , b"i) in Eq. (75)) are shown the white circles. The wavy lines denote the BFKL Pomerons. Fig. 3-c shows the Mueller diagram [51] for inclusive production of gluons.
in the following form [41]: where ρ t n and ρ p n are the parton densities in the target and projectile, respectively. These densities can be calculated from P n using Eq. (45). N BA is the scattering amplitude of two dipoles in the Born approximation of perturbative QCD (see Fig. 3). Eq. (75) simply states that we can consider the QCD parton cascade of Eq. (4) generated by the dipole of the size r for the c.m.f. rapidities from 0 to 1 2 Y , and the same cascade for the dipole of the size R, for the rapidities from 0 to − 1 2 Y . Generally speaking, for the dense system of partons at Y = 0 n-dipoles from upper cascade could interact with m dipoles from the low cascade, with the amplitude N m n {r i }, {r j } [41]. In Eq. (75) we assume that the system of dipoles that has been created at Y = 0 is not very dense. In this case and after integration over {r 1 } and {r j }, the scattering amplitude can be reduced to a system of enhanced BFKL Pomeron diagrams, which are shown in Fig. 3 The average number of dipoles at Y = 0 are determined by the inclusive cross section, which is given by the diagram of Fig. 3-c and which can be written at y → 0 as follows [52]: The average number of dipoles that enters the multiplicity distribution of Eq. (34) is equaln = N = 3 only if we assume that σ in ∼ Const. Indeed, the enhanced diagrams of Fig. 3-b lead to the inelastic cross section which is constant at high energy.

B. Hadron -hadron collisions
The first idea is to view a hadron as a dilute system of dipoles and use Eq. (77) for the average multiplicity, together with the multiplicity distribution of Eq. (34). However, the energy dependence of the mean multiplicity from Eq. (77) (n ∝ exp (∆ BFKL Y ) does not describe the experimental data (see Fig. 4). However, the experimental dependence of the mean multiplicity on energy can be parameterized asn ∝ exp (λ Y ), but with the value of λ = 0.1 ÷ 0.2 [54][55][56], which is far too small for ∆ BFKL . However, this power is close to the experimental behaviour of the deep inelastic structure function. Therefore, the main assumption of Ref. [5], that N ∼ xG(x, Q 2 ) , does not contradict the experimental data at least on the qualitative level [16].
On the other hand, the experimental data can be described in the framework of the CGC/saturation approach in which N BFKL were replaced by N BK [57]. Hence, we cannot view hadrons as the dilute system of dipoles, but rather have to consider them as the a dense system of dipoles. For such a situation we expect thatn ∝ Q 2 s (Y )/ᾱ S (see Refs. [24,57,[59][60][61]. Therefore the entanglement entropy in this case Note that for CGC approach of section III, the average multiplicity turns out to be proportional µ 2 ∼ Q 2 s , if the energy evolution is taken into account (see Ref. [24,30] for review). Frankly speaking, we do not have a theoretical tool to treat the dense-dense system scattering. For a general set of diagrams (see Fig. 5-a) we cannot use neither the Hamiltotian of Ref. [44], nor other theoretical methods. Therefore, our suggestion to use the multiplicity distribution of Eq. (34) withn determined by Eq. (77) with N BK is a conjecture. However, we can claim, that ρ 2 = ρ 1 (r 1 )ρ 1 (r 2 ) − ρ 1 (r 1 + r 2 ) (see Fig. 5-b and Fig. 5-c) , on the same theoretical grounds as the derivation of BK equation [27,28] since, the fact that the amplitude for two BFKL Pomeron production is equal to (N BK ) 2 , is used in the derivation of the BK equations.

C. Comparison with the experimental data
In order to compare the parton cascade with the experimental data, we first need to establish a relation between the multiplicity of hadrons with the multiplicities of partons, in the QCD parton cascade. Based on "parton liberation" picture [62] and on the " local parton-hadron duality" [63] we assume that S parton cascade = S hadrons . In other words, we suggest that there is no substantial entropy increase during the transformation of partons to hadrons. This relation is our hypothesis about confinement of quarks and gluons, and it has support in the fact, that the value of the entropy corresponds to a maximally entangled, equipartitioned state at a relatively modest average multiplicity of around N = 3 ÷ 6.
Hence, we use Eq. (34) for the hadron multiplicity distribution replacingn = N of partons, byn of hadrons. Following Ref. [5], we estimate, using Eq. (34), the value of the the cumulants where < · · · > denotes the average over the distribution in hadron multiplicity n. These quantities can be readily computed using Eq. (34)(see Ref. [5]). The result of these estimates are the following: C 2 = 2 − 1/n; C 3 = 6(n − 1)n + 1 n 2 ; Using the experimental multiplicity in the rapidity window |η| ≤ 0.5 equal[53] ton = 6.33 ± 0.46 at W=7 TeV and n = 3.72 ± 0.23 at W= 0.9 TeV we get from (Eq. (80)) the following predictions for the cumulants: C 2 1.83(1.73), C 3 5.08(4.46), C 4 18.6(15.31) and C 5 85.7(65.75). We put in parentheses the values at W = 0.9 T eV . The CMS experiments give (see Fig.6 Therefore, our estimates are in reasonably good agreement with the data, indicating that the parton distributions are close to the hadronic ones. Taking the limit ofn → ∞ we obtain the maximal values for the cumulants : C 2 = 2, C 3 = 6, C 4 = 24 and C 5 = 120 as a prediction for asymptotically high energies. Comparing these numbers to the experimental values listed above, we see that the multiplicity distribution measured at √ s = 7 TeV is already quite close to the expected asymptotic form.
In Fig. 6 and Fig. 7 we plot the multiplicity dependence in the form of the KNO scaling function [64]: One can see that at largen the distribution of Eq. (34) leads to and shows KNO scaling forn z. As far as we know, this is the first time that KNO scaling appears at ultra high energies on theoretical grounds for hadron-hadron collisions. At least, in the framework of the Pomeron calculus [65] KNO scaling is expected only for the intermediate range of energy [66,67].   In Fig. 6-a and Fig. 7-a we compare the CMS data [53] with the multiplicity distribution of Eq. (34). In spite of the good agreement at small z, one can see two major qualitative disagreements:(i) KNO scaling works better for |η| < 0.5 than at |η| < 2.4 in the data, while Eq. (34) predicts a different behaviour; and (ii) Eq. (34) predicts larger cross section in the region of large z, than is observed experimentally. The first disagreement is intimately related to the small value ofn at |η| < 0.5. In the framework of our approach the violation of KNO scaling will not be seen if we take N = 5.8 [68] at η = 0, instead of 6.33 = 30.4/4.8, which we used in these figures. The second disagreement is of a principle nature. As we have discussed we use Eq. (76),which is based on the assumption, that we do not have a dense system of parton at Y = 0. Certainly, such an assumption is not correct for events with large multiplicities. Therefore, this disagreement can be considered as an additional argument in the attempts to build a generalization of CGC approach to describe the hadron-hadron collisions.
In Fig. 6-b and Fig. 7-b we use the negative binomial distribution of Eq. (35) in which we fixed the parameter r from the moment n(n−1)

D. Back to QCD motivated parton model
In the previous section we inferred, that comparison with the experimental data indicates, that we cannot use Eq. (76), which stems from the assumption, that the system of partons, which is produced at Y = 0, is rather dilute. Unfortunately, we have not developed a theoretical approach in the framework of QCD to treat this problem. However, in the QCD parton model, that we have discussed above, a breakthrough has been made by [50], and an approach has been constructed, that sums all pomeron diagrams of the most general type (see Fig. 3 and Fig. 5 for examples).
In Ref. [50] a new Hamiltonian is suggested , which has the form where PM stands for " parton model".P and P are the BFKL Pomeron field in the model, where the sizes of the dipoles are fixed. This Hamiltonian in the limit of smallP reproduces the BK Hamiltonian ( see Ref. [50] and below for details). This condition is the most important one for fixing the form of H PM . The second of such conditions is that this Hamiltonian satisfies both t and s channel unitarity. γ in Eq. (83) has the physical meaning of the dipole-dipole scattering amplitude in the Born approximation of perturbative QCD and, being such, it is naturally small and of orderᾱ S . The most important ingredient of this approach is the generalization of the commutation relation, which has the form: Eq. (84) gives the correct factor (1 − γ)n that includes all multiple scattering corrections, while all the dipoles remain intact, and can subsequently scatter on an additional projectile or on target dipoles. For small γ, and in the regime where P andP are small themselves, we obtain consistent with our original expression. One can see that these commutation relations take into account the interaction of one dipole with many other partons and, therefore, we are going beyond the approximation, which is given by Eq. (76). Concluding this brief outline of this approach, we see, that for the first time we have a simple theory in which we can describe the interactions of dilute-dilute parton system scattering, as well as dilute-dense and dense-dense system interactions.
For H PM the cascade equation takes the form (see Eq. 5.8 of Ref. [50]): For small n (γ n < 1) one can see, that Eq. (86) reduces to Eq. (39). Hence, for such small n we have the multiplicity distribution of Eq. (41) with < n >= e ∆ Y . However, at large n Eq. (86) has the form We will show below that this equation gives the Poisson distribution with < n > = ∆ γ Y . Therefore, as we have guessed the interaction of one dipole with many dipoles at Y = 0 in Fig. 3 would lead to far fewer multiplicities than Eq. (41). Using Laplace transform of Eq. (27) and following the pattern, described in section II-C, we obtain the solution in ω-representation: where We have not found an elegant form for the inverse Laplace transform, but we can see the main qualitative features of this solution, assuming that for n < n 0 with γ n 0 ≈ 1 we have ω m = m ∆, but for n > n 0 ω m = ∆ γ . In this approach we obtain: In Fig. 8 we compare this multiplicity distribution with Eq. (34). One can see that at large multiplicities the modified distribution of Eq. (89a) and Eq. (89b) lead to the suppression of the parton emission, as we expected. Of course, this modified distribution is very approximate, and can only be used to clarify the qualitative features of the interaction of the partons in the exact approach.
To illustrate that the parton cascade of Eq. (86) is able to describe the experimental data. We calculate the first two P 1 and P 2 , taking integral over ω in Eq. (27): where ω 2 − ω 1 = ∆ − γ ∆ < ∆. P 2 in our notation with N = e ∆ Y can be re-written in the form: Assuming that the multiplicity distribution has the form: We can use this approximation, except for very large n. In Fig. 9 we compare Eq. (93) with the data. One can see that it provides quite a good description of the data.

V. CONCLUSSIONS
As has been discussed in the introduction, this paper has two main results. First, we prove that in QCD at high energies the multiplicity distribution has the form of Eq. (34), which was discussed in Ref. [5], in the framework of the parton model. We also show that the average number of gluons is not always related to the gluon structure function, and can depend both on energy and the size of the dipoles. However, the entanglement entropy is equal to S partons = ln N , where N is the average number of partons, confirming that the partonic state at high energy is maximally entangled [5]. In the case of DIS we prove that the average number of partons is related to the gluon structure function, only if we use the BFKL evolution equation for this structure function, can we prove that the multiplicity distribution has the form of Eq. (34).