\section{Introduction}
 
 
  During the past several years  a  robust 
relation
 between the principle features of high energy scattering and entanglement
 properties of the hadronic wave function   have been in focus of the high
 energy and nuclear physics communities  
\cite{KUT,PES,KOLU1,PESE,KHLE,BAKH,BFV,HHXY,KOV1,GOLE1,GOLE2,KOV2,NEWA,LIZA,FPV,TKU,KOV3,GOLE,KHLE1}.  In this paper, we continue to explore
 the relation between the entropy in the
 parton approach\cite{BJ,FEYN,BJP,Gribov} and  the entropy of
 entanglement in a proton wave function\cite{KHLE}.  In Ref.\cite{KHLE}, it is proposed that
 parton distributions can be defined in terms of the entropy of entanglement between the spatial region probed by deep inelastic scattering (DIS) and the rest of the proton.  This approach leads to a simple relation $S = \ln N $  between the average number of color-singlet dipoles and the entropy of the produced hadronic state $S$. This simple relation shows that  a proton becomes a maximally entangled state in the region of small Bjorken $x$.  All these conclusions were made from estimates in the simple, even naive  model for QCD cascade of color-singlet dipoles. However,  it has been demonstrated  in Refs.\cite{BAKH,GOLE1,GOLE2,TKU,KHLE1}  
 that these 
ideas are in qualitative and, partly, in  quantitative  agreement with the
 available experimental data.  Actually, it is shown in Ref.\cite{KHLE} that the simple cascade of color-singlet dipoles leads to the multiplicity distribution:

 
 \begin{equation} \label{MULTD}
\frac{\sigma_n}{\sigma_{ \rm in}}\,\,=\,\,\frac{1}{N}\,\left( \frac{N\,-\,1}{N}\right)^{n - 1}\,
\end{equation}
where $N$ is the average number of dipoles.

 
 

The goal of this paper is to study the multiplicity distribution and the
 entanglement entropy in the effective theory for  QCD at high energies
 (see Ref.\cite{KOLEB} for a general review). We have approached this problem in Refs.\cite{KHLE,GOLE} and have demonstrated that \eq{MULTD} arises in QCD cascades.  In    this paper we analyze the multiplicity distribution for Balitsky-Kovchegov (BK) cascade\cite{BK} in which one dipole at low energy generates a large number color-singlet dipoles at high energy. The equation for such a cascade  are known (see Refs.\cite{KOLEB,MUCD,LELU1,LELU2} ) and the first try to solve them have been undertaken in Ref.\cite{GOLE}. However, in this paper we return to this problem and study the multiplicity distribution using the new equation (Le, Mueller and Munier (LMM) equation)  for the probability generating function that has been derived in Ref.\cite{LMM}.
 
  In the next section we derive LMM equation from the equation for the BK parton cascade. In the rest of the paper we  discuss the equations for the factorial moments that follow from the LMM equation. We  show that every factorial moment satisfies the linear but inhomogeneous equation with the Balitsky, Fadin, Kuraev and Lipatov (BFKL) kernel\cite{BFKL,LI}. We  attempt  to solve these equations and  demonstrate that in the diffusion approximation to the BFKL kernel  factorial moments are equal to:
  
  \begin{equation} \label{C1}
M_k\,\,=\,\,k!\ N \,\Big( N \,\,-\,\,1\Big)^{k-1}
\end{equation}  
  
  We show that \eq{C1} leads to \eq{MULTD}.
  
  In section VI we suggest an approach  to go beyond    diffusion approximation, which cannot give a reliable description of the experimental data even in the leading order of perturbative QCD.  In this approach we propose  to solve exactly the    equations for the factorial moments and using the difference between the exact solution and  \eq{C1} ($\Delta M_k\,=\,M_k\left( \rm exact\right) \,-\,M_k\left( \eq{C1}\right)$) we develop  the way how to estimate the multiplicity distributions beyond diffusion approximation. 
   In conclusion section  we summarize our results.
 
 
 
 
 
 
 
 
 
  

\section{ General features of the  cascade of color-singlet dipoles in QCD }

In QCD at large number of colors $N_c$ ($N_c\,\,\gg\,\,1$) the color-singlet dipoles play the role of partons (see Ref.\cite{KOLEB} for review). 
 As discussed in Refs.\cite{KOLEB,MUCD,LELU1,LELU2}  for them we can write the following  equations:
  \begin{eqnarray}  \label{PC1}
&&\frac{\partial\,P_n\left(Y, \vec{r }, \vec{b};\,\vec{r}_1,\vec{ b}_1,\,\vec{r}_2 , \vec{b}_2\dots \vec{r}_i ,\vec{b}_i,
\dots \vec{r}_n, \vec{b}_n \right)}{ 
\partial\, Y }\,=\,-\,
\sum^n_{i=1}\,\omega_G(r_i) \,
P_n\left(Y, \vec{r }, \vec{b};\,\vec{r}_1,\vec{ b}_1,\,\vec{r}_2 , \vec{b}_2\dots \vec{r}_i ,\vec{b}_i,
\dots \vec{r}_n, \vec{b}_n \right) \\
&&~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~+\,\,\bar{\alpha}_S\,\sum^{n-1}_{i=1} \,\frac{(\vec{r}_i\,+\, 
\vec{r}_n)^2}{(2\,\pi)\,r^2_i\,r^2_n}\,
P_{n - 1}\left(Y, \vec{r},\vec{b};\,\vec{r}_1,\vec{b}_1,
\dots  (\vec{r}_i \,+\, \vec{r}_n), \vec{b}_{in},\dots \vec{r}_{n-1},\vec{b}_n
 \right)\nonumber
\end{eqnarray}
  where $P_n\left( Y ; \{r_i,b_i\}\right)$ is the probability to have $n$-dipoles
 of size $r_i$,  at impact parameter $b_i$ and  at rapidity $Y$\footnote
{ In the lab. frame rapidity $Y$ is equal to   $Y = y_{\rm dipole ~r} \,-\,y_{\rm dipoles~r_i}$, where $y_{\rm dipole ~r}$ is the rapidity of the incoming fast dipole and $y_{\rm dipole~r_i}$ is the rapidity of dipoles $r_i$. Note, that all rapidities of dipoles $r_i$ are the same in  \eq{PC1}.}
  . $\vec{b}_{in} $ in \eq{PC1} is equal to $\vec{b}_{in}
 \,=\,\vec{b}_i \,+\,\frac{1}{2} \vec{r}_i \,=\,\vec{b}_n \,-\,\frac{1}{2} \vec{r}_i$. $\omega_G\left( r\right)$ is defined below in  \eq{OMG}.
 
  
  \eq{PC1} 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
\begin{equation} \label{PCIC}
P_1 \left( Y  =  0,  \vec{r},\vec{b} ; \vec{r}_1,\vec{b}_1\right)\,\, =\,\,\,\delta^{(2)}\left( \vec{r}\,-\,\vec{r}_1\right)\,\delta^{(2)}\left( \vec{b}\,-\,\vec{b}_1\right);~~~~~~~P_{n>1}\left( Y  = 0; \{r_i\}\right) \,=\,0\end{equation}
which corresponds to the fact that we are discussing a dipole of 
 definite size which develops the parton cascade.
Since $P_n\left( Y ; \{r_i\}\right)$ is the probability to find dipoles $\{r_i\}$,
 we have the following sum rule 

\begin{equation} \label{SUMRU}
\sum_{n=1}^\infty\,\int \prod^n_{i=1} d^2 r_i \,d^2 b_i \,P_n\left( Y ; \{\vec{r}_i\,\vec{b}_i\}\right)\,\,=\,\,1 ,
\end{equation}
i.e. the sum of all probabilities is equal to 1.

 This QCD  cascade leads to Balitsky-Kovchegov (BK) equation 
\cite{BK,KOLEB} for the  amplitude and gives the theoretical
 description of the DIS.  We introduce the generating functional\cite{MUCD}

\begin{equation} \label{Z}
Z\left( Y, \vec{r},\vec{b}; [u_i]\right)\,\,=\,\,\sum^{\infty}_{n=1}\int P_n\left( Y,\vec{r},\vec{b};\{\vec{r}_i\,\vec{b}_i\}\right) \prod^{n}_{i=1} u\left( \vec{r}_i\,\vec{b}_i\right)\,d^2 r_i\,d^2 b_i
\end{equation}
 where $u\left( \vec{r}_i\,\vec{b}_i\right) \equiv\,u_i$ is an arbitrary function.
 The initial conditions of \eq{PCIC}  and the sum rules of \eq{SUMRU} take 
the following form for the functional $Z$:
\begin{subequations}
\begin{eqnarray}
Z\left( Y=0, \vec{r},\vec{b}; [u_i]\right) &\,\,=\,\,&u\left( \vec{r},\vec{b}\right);\label{ZIC}\\
Z\left( Y, r,[u_i=1]\right) &=& 1; \label{ZSR}
\end{eqnarray}
\end{subequations}

Multiplying both parts of \eq{PC1} by $\prod^{n}_{i=1} u\left( \vec{r}_i\,\vec{b}_i\right)$ and integrating over $r_i$ and $b_i$ we obtain the following linear functional equation\cite{LELU2};
\begin{subequations}
\begin{eqnarray}
&&\hspace{-0.7cm}\frac{\partial Z\left( Y, \vec{r},\vec{b}; [u_i]\right)}{\partial \,y} =\int d^2 r'\,  K\left( \vec{r}',\vec{r} - \vec{r'}|\vec{r}\right)\Bigg( - u\left( r, b\right)\,\,+\,\,u\left( \vec{r}',\vec{b} + \frac{1}{2}(\vec{r} - \vec{r}') \right) \,u\left( \vec{r} - \vec{r}',\vec{b} + \frac{1}{2}\vec{r}'\right)\Bigg) \frac{\delta\,Z}{\delta \,u\left( r, b \right)};\label{EQZ}\\
&&  K\left( \vec{r}',\vec{r} - \vec{r'}|\vec{r}\right)\,=\frac{1}{2 \,\pi}\frac{r^2}{r'^2\,(\vec{r} - \vec{r}')^2} ;\,~~~~~
 \omega_G\left( r\right)\,\,=\,\,\int d^2 r'  K\left( \vec{r}',\vec{r} - \vec{r'}|\vec{r}\right); \label{OMG}
 \end{eqnarray}
\end{subequations}
where $y\,\,=\,\,\bar{\alpha}_S\,Y$.

Searching for the solution of the form     $Z\left( [ u(r_i,b_i,Y)]\right)$  
for the initial conditions of \eq{ZIC}, \eq{EQZ} can be re-written as
 the non-linear equation \cite{MUCD}:
\begin{equation} \label{NEQZ}
\frac{\partial Z\left( Y, \vec{r},\vec{b}; [u_i]\right)}{\partial \,y}\,=\,\int d^2 r' K\left( \vec{r}',\vec{r} - \vec{r'}|\vec{r}\right)\Bigg\{Z\left( r',  \vec{b} + \frac{1}{2}(\vec{r} - \vec{r}'); [u_i]\right)  \,Z\left( \vec{r} - \vec{r'},  \vec{b} + \frac{1}{2}\vec{r}'; [u_i]\right)
\,\,-\,\,Z\left( Y, \vec{r},\vec{b}; [u_i]\right)\Bigg\}
\end{equation}
Therefore, the QCD parton cascade of \eq{PC1} takes into account 
 non-linear evolution. 
 
\section{Derivation of  Le, Mueller and Munier (LMM) equation }
In this section we derive the LMM equation which is  proposed in Ref.\cite{LMM}. First, we introduce the same notations as in Ref.\cite{LMM}:
 \begin{equation} \label{W}
 \tilde{w}_n\left( r, b, y\right)\,\,=\,\, \int \prod^n_{i=1} d^2 r_i \,d^2 b_i \,P_n\left( Y ; \{\vec{r}_i\,\vec{b}_i\}\right)
 \end{equation}
One can see that $ \tilde{w}_n\left( r, b, y\right)$ is the probability that the dipole  with size $r$ produces $n$ dipoles with all possible sizes. \eq{ZSR} reads as 
\begin{equation}\label{ZSRM}
 \sum_{n=1}^\infty \tilde{w}_n\left( r, b, y\right)\,\,=\,\,1
 \end{equation}


Taking all $u_i\left( r_i,b_i\right)\,\,=\,\,\lambda$ one can see that we can re-write  \eq{Z} in the form:
\begin{equation} \label{WLA} 
Z\left( Y, \vec{r},\vec{b}; [u_i = \lambda]\right)\,\,\equiv\,\tilde{w}_\lambda\left( \vec{r}, \vec{b}, y\right)\,\,=\,\, \sum_{n=1}^\infty  \lambda^n\,\,\tilde{w}_n\left( \vec{r}, \vec{b}, y\right)
\end{equation}

Plugging \eq{WLA} into \eq{NEQZ} we obtain the LMM equation in the form:
\begin{equation} \label{LMM}
\frac{\partial \tilde{w}_\lambda\left( \vec{r}, \vec{b}, y\right)}{\partial\,\,y}  =\int d^2 \,r' \, K\left( \vec{r}',\vec{r} - \vec{r'}|\vec{r}\right)\Bigg\{ \tilde{w}_\lambda\left( \vec{r}', \vec{b} + \frac{1}{2} \left( \vec{r}\,-\,\vec{r}'\right), y\right)\,\,\tilde{w}_\lambda\left( \vec{r}\,-\,\vec{r}', b\,+\,\frac{1}{2} \vec{r}', y\right)\,\,-\,\,\tilde{w}_\lambda\left( \vec{r},  \vec{b}, y\right)\Bigg\}
\end{equation}
  In addition, discussing multiplicity distribution $P_n \,=\,\frac{\sigma_n}{\sigma_{tot}}$, where $\sigma_n$ is the cross section for production of $n$  color-singlet dipoles, we need to integrate  $ w_n\left( r, b, y\right)$ over $b$. In this case the initial condition for the dipole cascade takes the form:
  \begin{equation} \label{PCIC1}
P_1 \left( Y  =  0,  \vec{r}; \vec{r}_1\right)\,\, =\,\,\,\delta^{(2)}\left( \vec{r}\,-\,\vec{r}_1\right)\;~~~~~~~P_{n>1}\left( Y  = 0; \{r_i\}\right) \,=\,0
\end{equation}  
 which leads to the probabilities, that do not depend on impact parameters.  Since in \eq{PC1} $b$  enters as a parameter,  $P_n\left(  Y; r_i\right)$, which does not depend on $b_i$, is also a solution to \eq{PC1},  which satisfies 
 \eq{PCIC1}. 
 
 
  \eq{LMM} reduces to
 \begin{equation} \label{LMM1}
\frac{\partial \tilde{w}_\lambda\left( \vec{r},  y\right)}{\partial\,\,y}  =\int d^2 \,r' \, K\left( \vec{r}',\vec{r} - \vec{r'}|\vec{r}\right)\Bigg\{ \tilde{w}_\lambda\left( \vec{r}', y\right)\,\,\tilde{w}_\lambda\left( \vec{r}\,-\,\vec{r}', y\right)\,\,-\,\,\tilde{w}_\lambda\left( \vec{r}, y\right)\Bigg\}
\end{equation} 
 
 \eq{LMM1} is a particular case of the general equation that has been derived in Ref.\cite{LMM}.  In this paper 
instead of \eq{W} the more general form of this equation is proposed, viz.:
\begin{equation} \label{WS}
 w_n\left( r, b,  y,  y_0 \right)\,\,=\,\, \int \prod^n_{i=1}\Big(  d^2 r_i \,d^2 b_i \,S\left( r_i, b_i, y_0\right)\Big) \,P_n\left( Y ; \{\vec{r}_i\,\vec{b}_i\}\right)
 \end{equation}
where $S\left( r_i, b_i, y_0\right)$ is the scattering S-matrix for elastic interaction of the dipole with size $r_i$  at rapidity $Y_0 $($y_0 \,=\,\bar{\alpha}_S\,Y_0$)  and at impact parameter $b_i$  with the target at $Y=0$. Since $S$ is a unitarity matrix, $ w_n\left( r, b, y; y_0\right)$ is the probability that the dipole  with size $r$ produces $n$ dipoles with all possible sizes, which interact with the target.
Bearing this in mind,  we see that \eq{ZSRM} holds for  $w_\lambda\left( r, b,  y; y_0\right)$, which is defined as

\begin{equation} \label{WLA} 
Z\left( Y, \vec{r},\vec{b}; [u_i = \lambda \,S\left( r_i, b_i, y_0\right)]\right)\,\,\equiv\,w_\lambda\left( \vec{r}, \vec{b}, y, y_0\right)\,\,=\,\, \sum_{n=1}^\infty  \lambda^n\,\,w_n\left( \vec{r}, \vec{b}, y, y_0\right)
\end{equation}


For the case of $P_n$ which do not depend on $b_i$, 
inserting in \eq{NEQZ} $u_i\left( r_i,b_i,\right) \,\,=\,\,\lambda S\left( r_i,b_i,y_0\right)$  we see that we obtain the  LMM equation in its original form (see Ref.\cite{LMM}):

\begin{equation} \label{LMMGEN}
\frac{\partial w_\lambda\left( \vec{r}, b,  y, y_0\right)}{\partial\,\,y}  = \int d^2 \,r' \, K\left( \vec{r}',\vec{r} - \vec{r'}|\vec{r}\right)\Bigg\{w_\lambda\left( \vec{r}', y, y_0\right)\,\,w_\lambda\left( \vec{r}\,-\,\vec{r}', y, y_0\right)\,\,-\,\,w_\lambda\left( \vec{r},  y, y_0\right)\Bigg\}
\end{equation}

\section{Average number of color-singlet  dipoles}
The average number of dipoles can be calculated using the following formula:
\begin{equation} \label{N}
N\left( \vec{r},  y, y_0 \right)\,\,=\,\,\langle| n |\rangle = \sum^\infty_{n=1} n\,w_n\left( r, y, y_0\right)\,=\,\frac{\partial\,\, w_\lambda\left( \vec{r},  y, y_0\right)}{\partial\,\,\lambda} \Big{|}_{\lambda\,=\,1}
\end{equation}

Differentiating \eq{LMM} with respect to $\lambda$ we obtain that
\begin{equation} \label{NEQN}
 \frac{\partial N\left( \vec{r}, y , y_0\right)}{\partial\,\,y}\,\,=\,\,\int d^2 r' K\left( \vec{r}',\vec{r} - \vec{r'}|\vec{r}\right)\Bigg\{N\left( \vec{r}', y, y_0 \right) \,+\,N\left( \vec{r}\,-\,\vec{r}', y, y_0 \right) \,\,-\,\,N\left( \vec{r}, y, y_0 \right)\Bigg\}
 \end{equation} 
 

 \eq{NEQN} shows that the average number of dipoles satisfies the linear  BFKL \cite{BFKL,LI} equation and 
 increases in the region of small $x$ (large $y$).  Therefore, we see that the general  QCD cascade  reproduces the main observation of Ref.\cite{KHLE} which was made in the oversimplified model  for the QCD cascade. In this model the dependence on the size of the dipoles were neglected. 
 
 The general solution takes the following form:

 \begin{equation} \label{SOLN}
N\left( \vec{r}, y, y_0 \right)\,\,\,=\,\,\,\int\limits^{\epsilon \,+\,i\,\infty}_{\epsilon \,-\,i\,\infty}\!\!\!\! \frac{d \gamma}{2\,\pi\,i}\,e^{ \chi\left( \gamma\right)\,y \,\,+\,\,\gamma\,\xi}\,n_{in}\left( \gamma, y_0\right)
\end{equation}
where  $\xi \,\,=\,\, \ln\left(\frac{1}{ r^2}\right)$  and $\chi(\gamma)$ is the BFKL kernel:

\begin{eqnarray} \label{BFKLKER}
\chi\left( \gamma \right)  =  2 \psi\left( 1\right) \,-\,\psi\left( \gamma\right)\,-\,\psi\left( 1 - \gamma\right)
= \left\{\begin{array}{l}\,\,\,\frac{1}{\gamma}\,\,\,\mbox{for}\,\,\,\gamma\,\to\,0  \leftarrow \,\,\mbox{double log approximation(DLA)};\\ \\
\underbrace{ 4 \ln 2}_{\omega_0} \,\,+\,\, \underbrace{14 \zeta\left( 3\right)}_{D} \left( \gamma - \frac{1}{2}\right)^2 \,\,\mbox{for}\,\,\,\gamma\,\to\,\frac{1}{2}\, \leftarrow \,\,\mbox{diffusion approximation(DA)};\\  \end{array}
\right.
\end{eqnarray}

 where $\psi(z)$ is the Euler $\psi$-function (see Ref.\cite{RY} formula{ \bf 8.36}). 
 
 $n_{in}\left( \gamma\right)$ has to be found from the initial condition $  N\left( \vec{r}, \vec{b}, y = 0,y_0=0 \right) \,\,=\,\,1$ (see \eq{PCIC} and \eq{ZIC}). It gives
 \begin{equation} \label{NIC}
 n_{in}\left( \gamma\right)\,\,=\,\,\frac{1}{\gamma}
 \end{equation}
 
 Introducing  multiplicity in the momentum representation:
 \begin{equation} \label{MMOM}
N\left( k_T, y\right) \,\,=\,\,\int d^2 r\,e^{- i \vec{k}_T \cdot \,\vec{r}} \frac{N\left( r,y\right)}{r^2},
\end{equation} 
 we can re-write \eq{NEQN}  in the form:
 
 \begin{equation} \label{NEQM}
\frac{\partial\,N\left( \vec{k}_T,  y,  y_0 \right)}{\partial\,y}\,\,=\,\, \int \frac{d^2 k_T}{(2 \pi)^2} K\left( \vec{k}_T,\vec{k}'_T\right)\,\,N\left( \vec{k}'_T, y, y_0 \right) 
\end{equation} 

where $K\left( \vec{k}_T,\vec{k}'_T\right)$ is the BFKL kernel in momentum representation:

\begin{equation} \label{KERM}
K\left( \vec{k}_T,\vec{k}'_T\right)\,\,N\left( \vec{k}'_T, y, y_0 \right)\,=\,\frac{1}{\left( \vec{k}_T - \vec{k}'_T\right)^2} \,\,N\left( \vec{k}'_T, y, y_0 \right)\,\,-\,\,\frac{k^2_T}{\left( \vec{k}_T - \vec{k}'_T\right)^2\,\left(\Lb \vec{k}_T - \vec{k}'_T\right)^2\,+\,k'^2_T\right) }\,\,N\left( \vec{k}_T, y, y_0 \right)
\end{equation}

 Solution to this equation has the same form of \eq{SOLN} but with the replacement of $\xi \,\to \xi' =    \ln k^2_T$ and
 \begin{equation} \label{NINM}
 n_{in}\left( \gamma\right)\,\,=\,\,\,\chi\left( \gamma\right)/\gamma
 \end{equation}
  \eq{NINM} reproduces  the value of $N\left( k_T, y=0, y_0=0\right)$, since \eq{MMOM} at $y=0$ leads to $ N\left( k_T, y\right) \,=\,\ln k^2_T + {\cal O}\left( 1/k^2_T\right)$. It is worth mentioning that  (1) in the double log approximation $ n_{in}\left( \gamma\right) = 1/\gamma^2$ , which leads to   $ N\left( k_T, y\right) \,=\,\ln k^2_T  $  and (2) in semi-classical approximation, which we will use below, we can neglect all corrections of the order of $1/k^2_T$.
 
 
\section{Equations for moments of the multiplicity distribution}
\subsection{The second moment}
\subsubsection{Equation}
We start  a derivation of the  evolution equation for the moments of the multiplicity distributions   considering the second moment, which has the following form:

\begin{equation} \label{M2}
M_2\left( \vec{r}, y,  y_0 \right)\,\,=\,\,\,\langle| n (n\,-\,1) |\rangle = \sum^\infty_{n=1} n\,(n\,-\,1) \,w_n\left( r, y, y_0\right)\,=\,\frac{\partial^2\,\, w_\lambda\left( \vec{r}, \ y, y_0\right)}{\partial\,\,\lambda^2} \Big{|}_{\lambda\,=\,1}
\end{equation}
 Taking the second derivative with respect to $\lambda$ from \eq{LMMGEN} we obtain the equation for $M_2\left( \vec{r},  y,  y_0 \right)$
 \begin{eqnarray} \label{M2EQ}
\frac{\partial\,M_2\left( \vec{r}, y,  y_0 \right)}{\partial\,y}\,\,&=&\,\, \int d^2 r' K\left( \vec{r}',\vec{r} - \vec{r'}|\vec{r}\right)\Bigg\{M_2\left( \vec{r}', y, y_0 \right) \,+\,M_2\left( \vec{r}\,-\,\vec{r}',  y, y_0 \right)\nonumber\\
 & \,\,+&\,\,2\,N\left( \vec{r}' , y, y_0 \right) \,N\left( \vec{r}\,-\,\vec{r}' , y ,y_0\right) -\,\,M_2\left( \vec{r},  y, y_0 \right)\Bigg\} 
\end{eqnarray}
 
 
 \eq{M2EQ} is a linear but inhomogeneous  equation with the  inhomogeneous  term, which is determined by the multiplicity of the dipoles.
 
 The initial conditions (see \eq{PCIC} and \eq{ZIC}) for this equation is
  \begin{equation} \label{M2IC}
  M_2\left( \vec{r},  y\,=\,y_0 , y_0 \right)  \,\,=\,\,0
  \end{equation}

First, lets us start to solve \eq{M2EQ} making first iteration at small $y = \Delta y$.

For $y=0$ $N\left( r, y =0\right) \,=\,1$ and $M_2\left( r, y =0\right)=0$, and hence 
\begin{equation} \label{M21I}
M^{(1)}_2 \,=\,2\,\Delta y \,\omega_G\left( r \right)
\end{equation}
 where $\omega  $ is given by \eq{OMG}. The first iteration of \eq{NEQN} leads to
$N^{(1)} \,=\,1\,\,+\,\,\,\Delta y \,\omega_G\left( r \right)$. Comparing these two estimates one can see that the first iteration can be written as the expansion of the solution $ M_2\left( r, y\right)\,\,=\,\,2\left( N^2\left( r, y\right)\,\,-\,\,
N\left( r, y \right)\Rb$ with respect to $\Delta \,y$. Hence, we see, that at small $\Delta y$ we obtain the simple expression for $M_2$, which turns out to be the same as for the multiplicity distribution of \eq{MULTD} for the simple toy model \cite{KHLE}.  We will try to prove this equation below, but we have  succeeded only in the  semi-classical  approximation.


\subsubsection{General solution}


The general solution to \eq{M2EQ} we can obtain going to momentum representation:

\begin{equation} \label{M2MOM}
m_2\left( k_T, y\right) \,\,=\,\,\int\frac{d^2 r}{(2 \pi)^2}e^{- i \vec{k}_T \cdot \,\vec{r}} \frac{M_2\left( r, y\right)}{r^2}
\end{equation}

\eq{M2EQ} takes the form:
\begin{equation} \label{M2EQMOM}
\frac{\partial\,m_2\left( \vec{k}_T,  y,  y_0 \right)}{\partial\,y}\,\,=\,\, \int \frac{d^2 k_T}{(2 \pi)^2} K\left( \vec{k}_T,\vec{k}'_T\right)\,\,m_2\left( \vec{k}'_T,  y, y_0 \right) 
  \,\,+\,\,2\,N^2\left( \vec{k}_T,  y, y_0 \right)  
\end{equation}

Taking the Mellin transform:
\begin{equation} \label{MLN}
   m_2\left( \vec{k}_T,  y,  y_0 \right) \,\,=\,\,\int\limits^{\epsilon + i  \infty}_{\epsilon - i \infty}\!\!\!\! \frac{d \gamma}{2\,\pi\,i}   e^{ \gamma\,\xi'}  \widetilde{m}_2\left( \gamma,y,y_0\right)~~~~~~\mbox{where}~~~\xi'\,\,=\,\,\ln k^2_T
   \end{equation} 
 from both part of \eq{M2EQMOM} we obtain:
 \begin{equation} \label{M2MLN}
 \frac{\partial\,\widetilde{m}_2\left( \gamma, y,  y_0 \right)}{\partial\,y}\,\,=\,\, \chi\left( \gamma\right) \widetilde{m}_2\left( \gamma,  y,  y_0 \right)\,\,+\,\,\,2\,\widetilde{n^2} 
\end{equation}
where $\chi\left( \gamma\right)$ is given by \eq{OMG} and $\widetilde{n^2} $ denotes the Mellin  image of $N^2\left( \vec{k}_T, \vec{b} , y, y_0 \right)$. 

\eq{M2MLN} has the following solution which satisfies the initial condition of \eq{M2IC}:

\begin{equation} \label{SOLMLN}
\widetilde{m}_2\left( \gamma,  y,  y_0 \right)\,\,=\,\, e^{\,\chi\left( \gamma\right)\,y} \int^y_0\!\!\!d y' \,2\,\widetilde{n^2}e^{\,-\,\chi\left( \gamma\right)\,y'}  +\,\,\widetilde{m}^{\rm BFKL}_2\left( \gamma,  y,  y_0 \right)
\end{equation}

where $\widetilde{m}^{\rm BFKL}_2\left( \gamma,  y,  y_0 \right)$ is a solution to the homogeneous  linear BFKL equation with the initial condition of \eq{M2IC}.  In the following, we neglect the contribution of  this term.~




\subsubsection{Semi-classical solution}


For large $y$ and $\xi$ we can use the semi-classical approximation( SCA,  see Refs.\cite{KOLEB,BKL} and references therein)  to take the integral over $y'$ in \eq{SOLMLN}. In this approximation we are searching for 
\begin{equation} \label{SCN}
N = e^{S_N} \,=\,e^{ \omega\left( \xi', y\right) \,y + \gamma\left( y,\xi'\right) \xi'}
\end{equation}
 where $\omega\left( \xi', y\right)$ and 
 $ \gamma\left( y,\xi'\right)$  are smooth functions of $y$ and $\xi'$:  \,$\partial \,  \omega'\left( \xi', y\right)/\partial\,y \,\ll\,\omega^2\left( \xi', y\right), \, \partial \, \omega\left( \xi', y\right)/\partial\,\xi'\,\ll\,\omega\left( \xi', y\right)\,\gamma\left( y,\xi'\right), \partial \,\gamma\left( \xi', y\right)/\partial\, \xi'\,\ll\,\gamma^2\left( \xi', y\right), \, \partial\, \gamma\left( \xi', y\right)/\partial \,\xi'\,\ll\,\omega\left( \xi', y\right)\,\gamma\left( y,\xi'\right)$. Such form  of $N$ stems from \eq{SOLN} if  we use the method of steepest descent for calculating the integral over $\gamma$.  Indeed, using this method one can see that 
 \begin{equation} \label{MOSD}
  \omega\left( \xi', y\right)\,\,=\,\,\chi\left(  \gamma_{\rm SP} \right) ;~~~~~\gamma\left( \xi', y\right) \,\,=\,\, \gamma_{\rm SP};  ~~~~\mbox{equation for $\gamma_{\rm SP}$:} ~~\frac{d \chi\left( \gamma\right)}{d \gamma}\Bigg{|}_{\gamma= \gamma_{\rm SP} }\!\!\!\!\!\!\!\!\!\!\!\!y\,\,=\,\,\xi';
  \end{equation}
  
  
  In the SCA the  Mellin image of $N^2$ can be written as follows:
  
  \begin{equation} \label{N2MLN}
    N^2\left( \xi',  y,  y_0 \right) \,\,=\,\,e^{2\,S_N} \,=\,e^{ 2\,\omega\left( \xi', y\right) \,y + 2\,\gamma\left( y,\xi'\right) \xi'}\,\,=\,\,   \int\limits^{\epsilon + i  \infty}_{\epsilon - i \infty}\!\!\!\! \frac{d \gamma}{2\,\pi\,i}   e^{\,\,2 \chi\left(\frac{1}{2} \gamma\right)\,y\,\,+\,\, \gamma\,\xi'}  H\left( \gamma\right)
   \end{equation} 
  Indeed, taking the integral by the method of steepest descent we obtain the following equation for the saddle point ($\gamma^{(2)}_{SP}$):
  \begin{equation} \label{MOSD1}
 2\frac{d\, \chi\left( \frac{1}{2}\gamma\right)}{d \gamma}\Bigg{|}_{\gamma= \gamma^{(2)}_{\rm SP} }y\,\,=\,\,\xi';  \end{equation}
 with the solution $ \gamma^{(2)}_{\rm SP}  \,=\,2\,\gamma_{SP} $, where $\gamma_{SP}$ is given by \eq{MOSD}. Plugging this solution into \eq{N2MLN} we see that 
  \begin{eqnarray} \label{N2MLN1}
    N^2\left( \xi',  y,  y_0 \right) \,\,&= &\,\, e^{\,\,2 \chi\left(\frac{1}{2} \gamma^{(2)}_{SP}\right)\,y\,\,+\,\, \gamma^{(2)}_{SP} \,\xi'}\int\limits^{\epsilon + i  \infty}_{\epsilon - i \infty}\!\!\!\! \frac{d \Delta \gamma}{2\,\pi\,i}   e^{\,\,\frac{d^2\chi\left( \frac{1}{2} \gamma^{(2)}_{SP} \right)}{d \gamma} \,y (\Delta \gamma)^2}  H\left(  \gamma^{(2)}_{SP} \right)\nonumber\\
    &=& e^{\,\,2 \chi\left( \gamma_{SP}\right)\,y\,\,+\,\,2\, \gamma_{SP} \,\xi'}\int\limits^{\epsilon + i  \infty}_{\epsilon - i \infty}\!\!\!\! \frac{d \Delta \gamma}{2\,\pi\,i}   e^{\,\,\frac{d^2\chi\left( \gamma_{SP} \right)}{d \gamma} \,y (\Delta \gamma)^2}  H\left(  2 \gamma_{SP} \right)   \end{eqnarray} 
 The integral over $\Delta \gamma$ leads to a smooth function, which in the SCA can be considered as a constant. Therefore, comparing \eq{N2MLN} and \eq{N2MLN1} one can see that \eq{N2MLN} is correct.  
 We can derive \eq{N2MLN} using a more general consideration.
  Actually, the  expression for the Mellin transform of $    N^2\left( \xi',  y,  y_0 \right) $ is the convolution in $\gamma$  of Mellin images of $N$, which has the following form:
\begin{equation} \label{MLNM21}
\widetilde{n^2}\,\,=\,\,\,\int\limits^{\epsilon + i\,\infty}_{\epsilon - i\,\infty}\!\! \frac{d \,\gamma'}{2\,\pi\,i}n_{in}\left( \gamma'\right)\,n_{in}\left( \gamma\,-\,\gamma'\right) e^{\left(  \chi\left( \gamma'\right) \,\,+\,\,\chi\left(\gamma\,-\, \gamma'\right)\Rb\,y}
\end{equation}
Taking the integral over $\gamma'$ using the method of steepest descent one can see that the equation for the saddle point has the following form:
\begin{equation} \label{SPGA} 
\frac{ d\,\chi\left( \gamma  - \gamma'\right)}{ d\,\gamma'} \,\,+\,\,\frac{ d\,\chi\left(  \gamma'\right)}{ d\,\gamma'}\,\,=\,\,0
\end{equation}

with the solution  $\gamma'_{SP} = \frac{1}{2} \gamma$. Plugging this solution in \eq{MLNM21} we reduce it to:
\begin{equation} \label{MLNM22}
\widetilde{n^2}\,\,=\,\,\,e^{ 2 \chi\left( \frac{1}{2} \gamma\right)\,y}\int\limits^{\epsilon + i\,\infty}_{\epsilon - i\,\infty}\!\! \frac{d \,\Delta \gamma'}{2\,\pi\,i}\,\,n^2_{in}\left( \frac{1}{2} \gamma\right) e^{\frac{d^2 \chi\left( \gamma'_{SP}\right)}{d \gamma'^2_{SP}}\,y\, \left( \Delta \gamma'\right)^2}
\end{equation}
which reproduces the  Mellin transform of \eq{N2MLN}.   
   
   
    Plugging in \eq{SOLMLN} $ \widetilde{n^2}\,\,=\,\,e^{ 2 \chi\left( \frac{1}{2} \gamma\right)\, y} H\left( \gamma\right)$ 
 we can take the integral over $y'$ and the solution has the form:
 
 \begin{equation} \label{SOLMLN1}
\tilde{m}_2\left( \gamma,  y,  y_0 \right)\,\,=\,\,\,e^{\,\,\chi\left( \gamma\right)\,y } H\left( \gamma\right) \int^y_0\!\!\!d y' \,2\,\widetilde{n^2}\,\,e^{\,-\,\chi\left( \gamma\right)\,y'}  
=\,\,\frac{2}{2  \chi\left( \frac{1}{2} \gamma\right)\,\,-\,\, \chi\left( \gamma\right)} \Bigg\{ e^{ 2 \chi\left(\frac{1}{2} \gamma\right)\, y} H\left( \gamma\right) \,\,-\,\, e^{ \chi\left( \gamma\right)\, y} H\left( \gamma\right)\Bigg\}
\end{equation}
Note, that in \eq{SOLMLN1} we neglected the contribution of  $\widetilde{m}^{\rm BFKL}_2\left( \gamma,  y,  y_0 \right)$ in \eq{SOLMLN}.

 Before fixing $H\left( \gamma\right)$
  we need to go back to coordinate representation. Indeed,  in this representation we have simple initial conditions for $M_2$  of \eq{M2IC}. Since all solutions are solutions of the linear equations  and $\gamma_{SP}\,\ll\,1$,  we can replace $\xi' =   \ln k^2_T$ by $\xi\,=\,\ln\frac{1}{ r^2}$.
   Bearing this in mind, we can reduce the solution to  the form:
  
 \begin{equation} \label{SOLMLN2}
M_2\left( \xi,  y,  y_0 \right)\,\,=\,\,\int\limits^{\epsilon + i\,\infty}_{\epsilon - i\,\infty}\!\!\!\frac{d\,\gamma}{2\,\pi\,i} \,\frac{2\,\chi\left( \gamma\right)}{2\,\chi\left( \frac{1}{2} \gamma\right)  \,\,-\,\, \chi\left( \gamma\right)} \Bigg\{ e^{ 2 \chi\left( \frac{1}{2} \gamma\right)\, y} \,\,-\,\, e^{ \chi\left( \gamma\right)\, y} \Bigg\}\,\frac{1}{\gamma} \,e^{ \gamma\,\,\xi}
\end{equation}
 
 First, we note, that taking the integral over $\gamma$ using the method of steepest descent, we reproduce \eq{SOLMLN1} with the particular choice of $H\left( \gamma\right)  =\chi\left( \gamma\right)/\gamma$, which  has been  discussed in \eq{NINM}.      Second, one can see that at $ y\,\to\,0$ this solution coincides with \eq{M21I}.
 
 In DLA  this solution takes the form (see \eq{BFKLKER}:
 \begin{subequations}  \begin{eqnarray} 
M^{\rm DLA}_2\left( \xi',  y,  y_0 \right)\,\,&=&\,\,\int\limits^{\epsilon + i\,\infty}_{\epsilon - i\,\infty}\!\! \frac{d \,\gamma}{2\,\pi\,i}\,\,\frac{2}{3} 
 \Bigg\{ e^{ \frac{4}{\gamma}\, y} \,\,-\,\, e^{ \frac{1}{\gamma}\, y} \Bigg\}\,\frac{1}{\gamma} \,e^{ \gamma\,\,\xi'}\label{SOLMLN31}\\
 &=& 
\frac{2}{3} \Bigg\{    N_{\rm SCA}^2\left( r,  y, y_0\right)\,\,-\,\,N_{\rm SCA}\left( r,  y, y_0\right) \Bigg\}\,\,\leftarrow \mbox{with semiclassical accuracy}\label{SOLMLN32}
\end{eqnarray}
 \end{subequations}  

 However, it turns out that \eq{SOLMLN31} reproduces at least two terms of the expansion at small values of 
 $y$ of the following relation:
 \begin{equation} \label{SOLMLN4} 
 M_2\left( \eq{SOLMLN31}\right)\,\,=\,\,2\Bigg( N^2\left( r,  y,y_0\right)\,\,-\,\,N\left( r,  y,y_0\right)\Bigg)\,\,=\,\,
 \,\,2\,y \xi \,\,+\,\,\frac{3}{2} \left( y\,\xi\right)^2\,+\,{\cal O}\left( \left( y\,\xi\right)^3\right) 
  \end{equation} 
  
  
  Concluding,  we see that  \eq{C1} does not hold in the DLA, but \eq{SOLMLN4} gives us some hope to find an approach in which it will be correct.
  
\subsubsection{Diffusion approximation}


Actually, the most adequate approach  at high energies is the diffusion one (see \eq{BFKLKER}).
Plugging in \eq{SOLMLN2} the BFKL kernel in the form $\chi\left( \gamma\right)\,\,=\,\,\omega_0 \,\,+D\,\left( \gamma\,-\,\frac{1}{2}\right)^2$  and taking the integral using the method of steepest  descent we obtain the saddle point 
\begin{equation} \label{SP}
\gamma_{SP} \,\,=\,\,\frac{1}{2} \,+ i \frac{\xi}{2\,D\,y}; ~~~ \,\,\,2\,\chi\left( \frac{1}{2} \gamma_{SP}\right)\,\,=\,\,2\,\omega_0\,\,-\frac{\xi^2}{2\,D\,y} ; ~~~~~\chi\left(  \gamma_{SP}\right)\,\,=\,\,\,\omega_0\,\,-\frac{\xi^2}{4\,D\,y} ;
\end{equation}
Considering $\frac{\xi}{2\,D\,y}\,\,\ll\,\,1$ we can neglect their contributions in the factor  $\frac{2\,\chi\left( \gamma\right)}{2\,\chi\left( \frac{1}{2} \gamma\right)  \,\,-\,\, \chi\left( \gamma\right)}$ reducing it to 2. Hence, \eq{SOLMLN2} reads 
\begin{equation} \label{SOLMLN6}
M^{\rm DA}_2\left( r,  y, y_0\right)\,\,=\,\,2\Bigg( N^2\left( r,  y, y_0\right)\,\,-\,\,N\left( r, y, y_0\right)\Bigg)
\end{equation}

We can obtain the solution of \eq{SOLMLN6} directly from the equation for $M_2$  (see \eq{M2EQ}), if we note that the BFKL kernel has   maxima at $r' \to 0$ and $|\vec{r} - \vec{r}'| \to 0$. In \fig{ker} we plot the term of \eq{M2EQ} , which is proportional to $N^2$:
\begin{equation} \label{I}
 \int d^2 r' K\left( \vec{r}',\vec{r} - \vec{r'}|\vec{r}\right)\,2\,N\left( \vec{r}' , y, y_0 \right) \,N\left( \vec{r}\,-\,\vec{r}' , y ,y_0\right)\,\, \propto\,\,  \int d^2 r'  K\left( \vec{r}',\vec{r} - \vec{r'}|\vec{r}\right)\,(\left( r'^2 \,\left( \vec{r} - \vec{r'}\right)^2\right)^{\frac{1}{2}} = \int d r' I\left( \tau =\frac{r'}{r} \right)
 \end{equation}
We can see from \fig{ker}, that $I(\tau)$ has a maximum at $\tau = 1$.  Note, that in   \eq{I} we introduce $\gamma = \frac{1}{2}$ , which corresponds to the DA, to estimates the  value of  this contribution.
 
   
     \begin{figure}[ht]
    \centering
  \leavevmode
      \includegraphics[width=10cm]{IvsK.pdf}  
      \caption{$I(\tau)$  of \eq{I} versus $\tau=r'/r$. }
\label{ker}
   \end{figure} 

 
 
 Bearing this observation in mind, we can re-write \eq{M2EQ} in the following form:
  \begin{eqnarray} \label{M2EQDA}
\frac{\partial\,M_2\left( \vec{r}, y,  y_0 \right)}{\partial\,y}\,\,&=&\,\, \int\limits_{r' < r} d^2 r' K\left( \vec{r}',\vec{r} - \vec{r'}|\vec{r}\right)  M_2\left( \vec{r}, y, y_0 \right) +  \int\limits_{r' < r} d^2 r'\, K\left( \vec{r}',\vec{r} - \vec{r'}|\vec{r}\right) M_2\left( \vec{r}', y, y_0 \right)\\
 &+& \int\limits_{|\vec{r} - \vec{r}'| < r}\!\!\!\!\! d^2 r'\, K\left( \vec{r}',\vec{r} - \vec{r'}|\vec{r}\right) M_2\left( \vec{r}\,-\,\vec{r}', y, y_0 \right) +  \int\limits_{|\vec{r} - \vec{r}'| < r}\!\!\!\!\!  d^2 r' \,K\left( \vec{r}',\vec{r} - \vec{r'}|\vec{r}\right) M_2\left( \vec{r}, y, y_0 \right)\nonumber\\
 & \,\,+&\,\,2\, \int\limits_{|\vec{r} - \vec{r}'| < r}\!\!\!\!\! d^2 r' \,K\left( \vec{r}',\vec{r} - \vec{r'}|\vec{r}\right) N\left( \vec{r} , y, y_0 \right) \,N\left( \vec{r}\,-\,\vec{r}' , y ,y_0\right)\nonumber\\
 & +&
 2 \int\limits_{r' < r} d^2 r'\, K\left( \vec{r}',\vec{r} - \vec{r'}|\vec{r}\right)  N\left( \vec{r} , y, y_0 \right) \,N\left( \,\vec{r}' , y ,y_0\right) 
   -\int d^2 r\, K\left( \vec{r}',\vec{r} - \vec{r'}|\vec{r}\right) \,\,M_2\left( \vec{r},  y, y_0 \right)\nonumber
\end{eqnarray}

One can see (1)  that in \eq{M2EQDA} the terms, which  are proportional to $M_2\left( \vec{r},  y, y_0 \right)$, cancel each other out, and (2) this equation can be presented in the form:
\begin{equation} \label{M2EQDA1}
\frac{\partial\,M_2\left( \vec{r}, y,  y_0 \right)}{\partial\,y}\,\,=\,\, \int d^2 r' K\left( \vec{r}',\vec{r} - \vec{r'}|\vec{r}\right)\Bigg\{M_2\left( \vec{r}', y, y_0 \right) \,\,+\,\,2\,N\left( \vec{r} , y, y_0 \right) \,N\left(\,\vec{r}' , y ,y_0\right)\Bigg\} \end{equation}


Substituting \eq{SOLMLN6} into \eq{M2EQDA1} we reduce it in the SCA to the equation:
\begin{equation} \label{M2EQDA2}
\frac{\partial\,2\left( N^2 - N\right)}{\partial\,y}\,=\, 2\,\chi\left( \gamma_{SP}\right) \left(  \,2 N^2  -  N\right)\,\,= \,\,2 \Big( \left( \chi \left( 2  \gamma_{SP}\right)\,+\,\chi\left( \gamma_{SP}\right) \right)\,N^2 - \chi\left( \gamma_{SP}\right)\,N\Big)
\end{equation}
Hence, one can see,  that in the limit of small $ \frac{\xi}{2\,D\,y}\,\,\ll\,\,1$,
$\chi \left( 2  \gamma_{SP}\right) = \chi \left(   \gamma_{SP}\right)$ and \eq{SOLMLN6} satisfies  \eq{M2EQDA1} and  \eq{M2EQDA2}.

  It is  worth noting that \eq{SOLMLN4} follows from the multiplicity distribution,  which is given by \eq{MULTD}.
We will concentrate our efforts on DA  in our presentation below.





\subsection{The third moment}

The third moment can be found from the generating function $\tilde{w}_\lambda\left( \vec{r}, \vec{b}, y, y_0\right)$ in the following way:

\begin{equation} \label{M3}
M_3\left( \vec{r},  y,  y_0 \right)\,\,=\,\,\,\langle| n (n\,-\,1)\,(n\,-\,2) |\rangle = \sum^\infty_{n=1} n\,(n-1)\,(n-2)w_n\left( r, y, y_0\right)\,=\,\frac{\partial^3\,\, w_\lambda\left( \vec{r},  y, y_0\right)}{\partial\,\,\lambda^3} \Big{|}_{\lambda\,=\,1}
\end{equation}

 Using \eq{M3}   we obtain the equation for $M_3\left( \vec{r},  y,  y_0 \right)$ from \eq{LMMGEN}:
 
 \begin{eqnarray} \label{M3EQ}
\frac{\partial\,M_3\left( \vec{r}, y,  y_0 \right)}{\partial\,y}\,\,&=&\,\, \int d^2 r' K\left( \vec{r}',\vec{r} - \vec{r'}|\vec{r}\right)\Bigg\{M_3\left( \vec{r}' , y, y_0 \right) \,+\,M_3\left( \vec{r}\,-\,\vec{r}', y, y_0 \right)\nonumber\\
 & \,\,+&\,\,6\,M_2\left( \vec{r}' , y, y_0 \right) \,N\left( \vec{r}\,-\,\vec{r}', y ,y_0\right) -\,\,M_3\left( \vec{r}, y, y_0 \right)\Bigg\} \nonumber
\end{eqnarray}
Rewriting \eq{M3EQ} in  momentum representation(see \eq{M2MOM})  we reduce  this eqaution to  the form:
\begin{equation} \label{M3EQMOM}
\frac{\partial\,m_3\left( \vec{k}_T,  y,  y_0 \right)}{\partial\,y}\,\,=\,\, \int \frac{d^2 k'_T}{(2 \pi)^2} K\left( \vec{k}_T,\vec{k}'_T\right)\,\,m_3\left( \vec{k}'_T,  y, y_0 \right) 
  \,\,+\,\,6\,m_2\left( \vec{k}_T, y, y_0 \right)\, n\left( \vec{k}_T, y, y_0 \right)   
\end{equation}
where we use lowcase letters denoting the moments in the momentum representation.

The Mellin image of \eq{M3EQMOM}  has the form:
 \begin{equation} \label{M3MLN}
 \frac{\partial\,\widetilde{m}_3\left( \gamma,  y,  y_0 \right)}{\partial\,y}\,\,=\,\, \chi\left( \gamma\right) \widetilde{m}_3\left( \gamma,  y,  y_0 \right)\,\,+\,\,\,6 \,\widetilde{n^2\,n} 
\end{equation}
where $\widetilde{n^2\,n} $ is the Mellin image of $m_2\left( \vec{k}_T,  y, y_0 \right)\, N\left( \vec{k}_T , y, y_0 \right)$ which has the general form:

 \begin{equation} \label{M2MLN70}
\widetilde{n^2\,n}\left( \gamma, y\right)\,\,=\,\,\,\int\limits^{\epsilon + i\,\infty}_{\epsilon - i\,\infty} \frac{d \,\gamma'}{2\,\pi\,i}\,
\widetilde{n^2}\left( \gamma', y\right)\,\widetilde{n}\left( \gamma  -  \gamma', y\right)\end{equation}
where $\widetilde{n^2}\left( \gamma, y\right)$ is determined by (see \eq{SOLMLN2})
\begin{equation} \label{M2MLN7}
\widetilde{n^2}\,\,=\,\,\frac{2\,\chi\left( \gamma\right)}{2\,\chi\left( \frac{1}{2} \gamma\right)  \,\,-\,\, \chi\left( \gamma\right)} \Bigg\{ e^{ 2 \chi\left( \frac{1}{2} \gamma\right)\, y} \,\,-\,\, e^{ \chi\left( \gamma\right)\, y} \Bigg\}\,\frac{1}{\gamma}
\end{equation}
Plugging \eq{SOLN} and \eq{M2MLN7} into \eq{M2MLN70}   one can see that $\widetilde{n^2\,n}\left( \gamma, y\right)$ is a sun of two terms.  
Taking the integral over $\gamma'$ in each of them in the saddle point approximation  (see \eq{MLNM21} - \eq{MLNM22}) we obtain  two equations for the saddle points:
\begin{equation} \label{M2MLN71}
2\frac{d\,\chi\left( \frac{1}{2} \gamma'\right)}{d\gamma' } + \frac{d\,\chi\left( \gamma -  \gamma'\right)}{d\gamma' } = 0~~~\mbox{with} ~~~\gamma'_{SP} = \frac{2}{3} \gamma;  ~~~~~
\frac{ d\,\chi\left( \gamma  - \gamma'\right)}{ d\,\gamma'} \,\,+\,\,\frac{ d\,\chi\left(  \gamma'\right)}{ d\,\gamma'}\,\,=\,\,0;~ ~~\mbox{with}\,\gamma'_{SP} = \frac{1}{2}  \gamma;
\end{equation}

Hence,  $\widetilde{n^2\,n}\left( \gamma, y\right)$  takes the form:
 \begin{equation} \label{MLNM211}
\widetilde{n^2\,n}\left( \gamma, y\right)\,\,=\,\,\,\frac{2\,\chi\left( \gamma\right)}{2\,\chi\left( \frac{1}{3} \gamma\right)  \,\,-\,\, \chi\left(\frac{1}{2} \gamma\right)} \Bigg\{ e^{ 3\chi\left( \frac{1}{3} \gamma\right)\, y} \,\,-\,\, e^{ \chi\left( \frac{1}{2} \gamma\right)\, y} \Bigg\}\,H_2\left( \gamma\right)
\end{equation}
The solution does not depend on the  form of $H_2\left( \gamma\right)$ , which we will specify below.

Plugging this equation in \eq{M3MLN} we obtain the solution:
\begin{eqnarray} \label{SOLM3MLN}
\widetilde{m}_3\left( \gamma,  y,  y_0 \right)\,\,&=&\,\,6\left(\frac{1}{3 \,\chi\left( \frac{1}{3} \gamma\right) \,\,-\,\,\chi\left(  \gamma\right)}\right)\,\left(\frac{2\,\chi\left( \gamma\right)}{2\,\chi\left( \frac{1}{3} \gamma\right)  \,\,-\,\, \chi\left(\frac{1}{2} \gamma\right)}\right) \Bigg\{ e^{ 3\,\chi\left( \frac{1}{3} \gamma\right)\,y} \,\,-\,\, e^{ \,\chi\left( \gamma\right)\,y}\Bigg\}\,\,\nonumber\\
 &  -&\,\,\left( \frac{1}{2\,\chi\left( \frac{1}{3} \gamma\right)  \,\,-\,\, \chi\left(\frac{1}{2} \gamma\right)}\right)\left( \frac{2\,\chi\left( \gamma\right)}{2\,\chi\left( \frac{1}{2} \gamma\right)  \,\,-\,\, \chi\left( \gamma\right)}\right)\Bigg\{ e^{ 2\,\chi\left( \frac{1}{2} \gamma\right)\,y} \,\,-\,\, e^{ \,\chi\left( \gamma\right)\,y}\Bigg\}\,H_2\left( \gamma\right)
 \end{eqnarray}
 Going to the coordinate representation as was discussed above  and choosing $H_2\left( \gamma\right)\,=\,\chi\left( \gamma\right)/\gamma$, which reproduces the correct initial conditions we obtain 
  \begin{eqnarray} \label{SOLM3MLN2}
M_3\left( \xi, y,  y_0 \right)\,\,&=&\,\,6\,\int\limits^{\epsilon + i\,\infty}_{\epsilon - i\,\infty}\frac{d\,\gamma}{2\,\pi\,i} ,\left(\frac{\chi\left( \gamma\right)}{3 \,\chi\left( \frac{1}{3} \gamma\right) \,\,-\,\,\chi\left(  \gamma\right)}\right)\,\left(\frac{2\,\chi\left( \gamma\right)}{2\,\chi\left( \frac{1}{3} \gamma\right)  \,\,-\,\, \chi\left(\frac{1}{2} \gamma\right)}\right) \Bigg\{ e^{ 3\,\chi\left( \frac{1}{3} \gamma\right)\,y} \,\,-\,\, e^{ \,\chi\left( \gamma\right)\,y}\Bigg\}\,\,\nonumber\\
 &  -&\,\,\left( \frac{\chi\left( \gamma\right)}{2\,\chi\left( \frac{1}{3} \gamma\right)  \,\,-\,\, \chi\left(\frac{1}{2} \gamma\right)}\right)\left( \frac{2\,\chi\left( \gamma\right)}{2\,\chi\left( \frac{1}{2} \gamma\right)  \,\,-\,\, \chi\left( \gamma\right)}\right)\Bigg\{ e^{ 2\,\chi\left( \frac{1}{2} \gamma\right)\,y} \,\,-\,\, e^{ \,\chi\left( \gamma\right)\,y}\Bigg\}\,\frac{1}{\gamma}e^{ \gamma\,\,\xi}
\end{eqnarray} 

Using the method of steepest descent and neglecting contributions of the order of  $\frac{\xi'}{2 \,D\,y} $  in all  pre-exponential factors,  we see that
\begin{equation} \label{SOLM3MLN3}
M_3\left( \xi,  y,  y_0 \right)\,\,=\,\,6 \,N\left( \xi, y,  y_0 \right)\,\left( N\left( \xi,  y,  y_0 \right)\,\,-\,\,1\right)^2,
\end{equation}
which is the same as for the multiplicity distribution of \eq{MULTD}.


\subsection{General approach}
The equation for a general  moment  
\begin{equation} \label{MK}
M_k\left( r,  y, y_0\right) = \langle| n (n\,-\,1) \dots (n - k+1) |\rangle = \sum^\infty_{n=1}\left(  n\,(n\,-\,1) \dots (n-k+1)\right)\,w_n\left( r,  y, y_0\right)\,=\,\frac{\partial^k\,\, w_\lambda\left( \vec{r},  y, y_0\right)}{\partial\,\,\lambda^k} \Big{|}_{\lambda\,=\,1}
\end{equation}
we can obtain by  differentiating \eq{LMMGEN} with respect to $\lambda$.  It has the simple form in the momentum representation (see \eq{M2MOM}):
\begin{equation} \label{MKEQ}
\frac{\partial}{\partial\,y}m_k\left(  k_T, y, y_0\right)\,\,=\,\,\int \frac{d^2 k_T}{(2 \pi)^2} K\left( \vec{k}_T,\vec{k}'_T\right)\,\,m_k\left( \vec{k}'_T,  y, y_0 \right) 
  \,\,+\,\,\sum^{k-1}_{i=1}\,\frac{k!}{(k - i)!\,i!}m_{k - i}\left( \vec{k}_T, \ y, y_0 \right)\, m_{i}\left( \vec{k}_T, y, y_0 \right)  
\end{equation}

As we have discussed, we can go back to coordinate representation  in \eq{MKEQ}, since $m_k$ are solution to the linear equations, and, therefore,  have the Mellin image 

\begin{equation} \label{MKMLN}
 m_k\left( \vec{k}_T, y \right)\,\,\,=\,\,\,\int\limits^{\epsilon \,+\,i\,\infty}_{\epsilon \,-\,i\,\infty}\!\! \frac{d \gamma}{2\,\pi\,i}\,e^{ \omega\left( \gamma\right)\,y \,\,+\,\,\gamma\,\xi'}\,m_{k0}\left( \gamma\right)~~~~\mbox{with}~~~\xi' \,\,=\,\,\ln k^2_T
 \end{equation}
In coordinate representation we have
\begin{equation} \label{MKMLN1}
 M_k\left( r,  y \right)\,\,\,=\,\,\,\int\limits^{\epsilon \,+\,i\,\infty}_{\epsilon \,-\,i\,\infty}\!\! \frac{d \gamma}{2\,\pi\,i}\,e^{ \omega\left( \gamma\right)\,y \,\,+\,\,\gamma\,\xi}\,m_{k0} \left( \gamma\right)\,H_k(\gamma)~~~~\mbox{with}~~~\xi \,\,=\,\,\ln\left( \frac{1}{ r^2}\right)
 \end{equation}
The coordinate image of 
$m_{k - i}\left( \vec{k}_T, \ y, y_0 \right)\, m_{k}\left( \vec{k}_T, y, y_0 \right)$ is  $
\int d^2 r'\,K\left( \vec{r}',\vec{r} - \vec{r'}|\vec{r}\right)\,M_{k-i}\left( \vec{r}',  y, y_0 \right) \,M_i\left( \vec{r}\,-\,\vec{r}',y, y_0\right)$. On the other hand, in SCA   $m_{k - i}\left( \vec{k}_T, \ y, y_0 \right)\, m_{i}\left( \vec{k}_T, y, y_0 \right) $ has the image in $\gamma$-representation, which is equal to ${\rm Const} \exp\left( k \chi\left( \frac{1}{k}\gamma\right) y\right)\,m_{k-i,0}\left( \frac{\gamma}{k}\right)\,m_{i,0}\left( \frac{\gamma}{k}\right)\,\widetilde{H}_k\left( \gamma\right)$   (see \eq{SOLMLN1} and \eq{MLNM211})\footnote{ Actually, the $\gamma$ image of $m_k$
is a sum of the terms with different $k$ (see \eq{SOLMLN1} and \eq{MLNM211}), but  this does not change  the conclusion.}. Using this image,  one can see, that the coordinate representation for $m_{k - i}\left( \vec{k}_T, \ y, y_0 \right)\, m_{k}\left( \vec{k}_T, y, y_0 \right)$ can be reduced to $\int d^2 r' K\left( \vec{r}',\vec{r} - \vec{r'}|\vec{r}\right)\,\,M_{k - i}\left( r', y, y_0 \right)\, M_{i}\left( r',  y, y_0 \right)$, which means, that $\widetilde{H}_k\left( \gamma\right)\,=\,\chi\left( \gamma\right)$ as we expect from \eq{SOLMLN2} and \eq{SOLM3MLN2}.

Finally we can re-write
\eq{MKEQ} in the coordinate representation in the form:
\begin{eqnarray} \label{MKEQ2}
\frac{\partial}{\partial\,y} M_k\left(  r, y, y_0\right)\,\,&=&\,\, \int d^2 r' K\left( \vec{r}',\vec{r} - \vec{r'}|\vec{r}\right)\,
\Bigg\{M_k\left( \vec{r}',  y, y_0 \right) \,+\,M_k\left( \vec{r}\,-\,\vec{r}', y, y_0 \right)\,-\,M_k\left( \vec{r} ,y, y_0 \right) \nonumber\\
& +&\,\,\,\sum^{k-1}_{i=1}\,\frac{k!}{(k-i)!\,i!}M_{k - i}\left( r', y, y_0 \right)\, M_{i}\left( r',  y, y_0 \right)\Bigg\} 
\end{eqnarray}
  .
                                                                                                                                            

Assuming that for all $i\,\leq\,k - 1$ 
\begin{equation} \label{MI}
M_i\left( r , y, y_0 \right)\,\,=\,\,i! \,N\left( r, y, y_0 \right)\,\left( N\left( r, y, y_0 \right)\,-\,1\right)^{i - 1},
\end{equation} 
 which follows from \eq{MULTD}, we will prove that for $i=k$ we have the same expression.


Plugging \eq{MI} in  \eq{MKEQ2} we get the inhomogeneous term in  the form $k! \left( k - 1\right) N^2\left( N \,-\,1\right)^{k-2}$. In the following we will use that the Mellin  image of $N^i\left( r, y, y_0 \right)$ ($\widetilde{n^i}$) is  equal to
\begin{equation} \label{MIMLNNI}
\widetilde{n^i}\,\,=\,\,e^{i\,\chi\left( \frac{1}{i}\gamma\right)\,y} \frac{1}{\gamma}\,\,=\,\,\,\int\limits^{\epsilon \,+\,i\,\infty}_{\epsilon \,-\,i\,\infty} \frac{d \omega}{2\,\pi\,i}\,e^{\omega \,y} \frac{1}{\omega\,\,-\,\,i\,\chi\left( \frac{1}{i}\gamma\right)}\frac{1}{\gamma};~~~~~~~\mbox{double Mellin image} \,\,\widetilde{\widetilde{n^i}}\,\,=\,\,\frac{1}{\omega\,\,-\,\,i\,\chi\left( \frac{1}{i}\gamma\right)}\frac{1}{\gamma};\end{equation}
which can be derived using the method of steepest descent  in  the estimates of the integrals over $\gamma$'s (see \eq{MLNM21} for example).

In the double Mellin transform \eq{MKEQ2} takes the form:
\begin{equation} \label{MIEQ2MLN}
\left( \omega\,\,-\,\,\chi\left( \gamma\right)\Rb\widetilde{\widetilde{M_k}}\,\,=\,\,k!\,(k - 1)\chi\left( \gamma\right)\sum^{k-2}_{l=0} 
\frac{(-1)^l (k - 2)!}{(k - 2 - l)! \,l!}\,\frac{1}{\omega\,-\,\left( k - l\right) \chi\left(\frac{1}{k - l}\gamma\right)}\frac{1}{\gamma}
\end{equation}

Hence from \eq{MIEQ2MLN} we have:
\begin{eqnarray} \label{SOLMK}
&&M_k\left(  r,y, y_0\right)\,\,=\,\,\,k!\,(k - 1)\sum^{k-2}_{l=0} \frac{(-1)^l (k - 2)!}{(k - 2 - l)! \,l!}\,\\
&& \times\,\,\int\limits^{\epsilon \,+\,i\,\infty}_{\epsilon \,-\,i\,\infty} \frac{d \omega}{2\,\pi\,i}\,\int\limits^{\epsilon \,+\,i\,\infty}_{\epsilon \,-\,i\,\infty} \frac{d \gamma}{2\,\pi\,i}\,e^{\omega \,y\,\,+\,\,\gamma\,\xi'} 
\frac{\chi\left( \gamma\right)}{(n - l) \chi\left(\frac{1}{k - l}\gamma\right)\,-\,\chi\left( \gamma\right)}\Bigg\{
\frac{1}{\omega\,-\,\left( k - l\right) \chi\left(\frac{1}{k - l}\gamma\right)}\,\,-\,\,\,\frac{1}{\omega\,-\,\chi\left( \gamma\right)}\Bigg\}\frac{1}{\gamma}\nonumber\\
&&=\,k!\,(k - 1)\sum^{n-2}_{l=0} \frac{(-1)^l (k - 2)!}{(k - 2 - l)! \,l!}\int\limits^{\epsilon \,+\,i\,\infty}_{\epsilon \,-\,i\,\infty} \frac{d \gamma}{2\,\pi\,i}\,e^{\,\,\gamma\,\xi'} 
\frac{\chi\left( \gamma\right)}{(k - l) \chi\left(\frac{1}{k - l}\gamma\right)\,-\,\chi\left( \gamma\right)}\Bigg\{
e^{ \left(  k - l\right) \chi\left(\frac{1}{k - 1}\gamma\right)\,y}\,\,-\,\,\,e^{\chi\left( \gamma\right)\,y}\Bigg\}\frac{1}{\gamma}\nonumber\end{eqnarray}
In \eq{SOLMK} we  calculate  the integrals over $\omega$ closing the contour of integration on  poles.
Taking the integral over $\gamma$ using the method of steepest descent  in the diffusion approximation  and neglecting the corrections of the order of $\xi'/Dy$ (see discussions above)  we reduce \eq{SOLMK}  to the following expression:
\begin{eqnarray} \label{SOLMK1}
M_k\left(  r,  y, y_0\right)\,\,&=&\,\,k!\,\sum^{k-2}_{l=0} \frac{(-1)^l (k - 1)!}{(k - 1 - l)!\, l!}\left( N^{k-l}\left(  r,  y, y_0\right)\,\,-\,\,N\left(  r, b, y, y_0\right)\Rb\,\,\nonumber\\
&=& k!\,\sum^{k-1}_{l=0} \frac{(-1)^l (k - 1)!}{(k - 1 - l)! \,l!} N^{k-l}\left(  r,  y, y_0\right) \,=\,   k!\, N\left(  r,  y, y_0\right)\,\left( N\left(  r,  y, y_0\right)\,\,-\,\,1\right)^{k - 1}
 \end{eqnarray}

Since we have obtained \eq{MI} for $i=2$ and $i=3$ , we prove this equation for any value of $i$.




\section{Finding corrections and comparison with experiments}

\eq{MI} generates the multiplicity distribution of    \eq{MULTD}. However, several questions arise,  when we wish to compare this distribution with the experimental data, let say with DIS. The average multiplicity of the color-singlet dipoles is equal to the sea quark structure function $x\Sigma_{\rm sea}\left( x, Q^2\right)$\cite{KHLE1}. On the other hand, in \eq{MI} the multiplicity enters in the coordinate representation.  In diffraction approximation the momentum and coordinate representations are related by the replacement $\ln k^2_T\,\to -\ln r^2$. Therefore, my suggestion is to use \eq{MULTD} with $N\,=x\Sigma_{\rm sea}\left( x, Q^2\right)$ but to calculate the corrections to this distribution.  The multiplicity  distribution can generally be written, using the cumulant generating function $f\left( \lambda\right)$  as follows \cite{MUMULT,MGF}:
\begin{equation}
\frac{\sigma_n}{\sigma_{ \rm in}}\,\,=\,\,\oint\frac{d \lambda}{2\,\pi\,i} \frac{e^{ f\left( \lambda\right)}}{\lambda^{n + 1}}
\end{equation}
where the contour of  integration  is the circle around the point  $\lambda =0$  and  $f\left( \lambda\right)$ is the cumulant generating function, which  is defined as
\begin{equation} \label{CGF}
f\left( \lambda\right)\,\,=\,\,\sum^\infty_{n=1} \frac{\kappa_n }{n!} \left( \lambda - 1\right)^n
\end{equation}
where $\kappa_n$ are cumulants. Generally speaking, we have the following definition for the cumulants:
\begin{equation} \label{CUM}
\kappa_1\,=\,N;~~~\kappa_2\,=\,M_2\,-\,N^2;~~~ \kappa_3\,=\,M_3\,-\,3\,M_2 \,N\,+\,2\,N^2;~~~\kappa_4\,=\,M_4\,-\,4\,M_3\,N\,-\,3\,M^2_2\,+\,12\,M_2\,N^2\,-\,6\,N^4;
\end{equation}
where $M_k$ are the factorial moments, that we have discussed above (see \eq{MI}).

In our case we  can view $f\left( \lambda\right)$ as a sum of $f\left( \lambda\right)\,\,=\,\,f^{\eq{MULTD}}\left( \lambda\right)\,\,+\,\,\Delta f\left( \lambda\right)$. $f^{\eq{MULTD}}\left( \lambda\right)$ generates the multiplicity distributions of \eq{MULTD}, which includes the most dominant contributions and ,in particular, the average number of color-singlet dipoles is taken into account exactly.  We suggest to introduce function 
$\Delta f\left( \lambda\right)$ in the following way:
\begin{equation} \label{CGFD}
\Delta f\left( \lambda\right)\,\,=\,\,\sum^\infty_{n=1} \frac{\Delta\kappa_n }{n!} \left( \lambda - 1\right)^n
\end{equation}
with
\begin{equation} \label{DCUM}
\Delta \kappa_1\,=\,0;~~\Delta \kappa_2\,\,=\,\,M_2\left( k_T, y\right) \,-\,2\,N\,\left( N\,- \,1\right) ;~~~~
 \Delta \kappa_3\,=\,M_3\left( k_T, y\right)\,-  6\,N \left(  N\,-\,1\right)^2 \, -\,3\,\Delta \kappa_2 \,N\,;\end{equation}
where $M_2\left( k_T, y\right)$  and $M_3\left( k_T, y\right)$ are  the exact solutions to \eq{M2EQMOM} and  \eq{M3EQMOM}, respectively.


Introducing two multiplicity distributions:
\begin{equation} \label{2D}
P_n\left( N\right)\,\,=\,\,\,\,\oint\frac{d \lambda}{2\,\pi\,i} \frac{e^{ f^{\eq{MULTD}}\left( \lambda\right)}}{\lambda^{n + 1}};~~~~~~~\widetilde{P}_n\left( N\right)\,\,=\,\,\,\,\oint\frac{d \lambda}{2\,\pi\,i} \frac{e^{\Delta f\left( \lambda\right)}}{\lambda^{n + 1}};
\end{equation}
one can see that the resulting multiplicity distribution takes the form

\begin{equation} \label{RMULTD}
\frac{\sigma_n}{\sigma_{ \rm in}}\,\,=\,\,\sum^n_{k=0}\,\frac{n!}{(n - k)!\,k!}\,P_{n - k}\left( N\right)\, \widetilde{P}_k\left( N\right)
\end{equation}

In the case of $\Delta \kappa_2\,\neq\,0$ but  $\Delta \kappa_n\,=\,0$ for $n\,>\,2$ the distribution $\widetilde{P}_n\left( N\right)$ has been found in Ref.\cite{MUMULT} and it has the following form in our notations:

\begin{equation} \label{TILDEP}
\widetilde{P}_n\left( \Delta \kappa_2\right) \,\,=\,\,e^{ \frac{\Delta \kappa_2}{2}}\frac{\left( - i \sqrt{ \frac{\Delta \kappa_2}{2}}\right)^n}{n!}\,{\rm H_n}\left( i \sqrt{ \frac{\Delta \kappa_2}{2}}\right)
\end{equation}
where ${\rm H}_n $ is   Hermite polynomial (see formula {\bf 8.95} in Ref.\cite{RY}).


\section{Conclusions}
 
This paper has two main 
results. 
First, we derived the BFKL linear, inhomogeneous equation for the factorial moments of  multiplicity distribution ($M_k$) from LMM  equation. In particular, the equation for the average multiplicity of the color-singlet ($N$)
turns out to be the homogeneous  BFKL equation which leads to the power-like growth in the region of small $x$.  From these equations it follows that $M_k\,\,\propto\,\,N^k$ at small $x$. 


Second, using the diffusion approximation for the BFKL kernel, which is generally  considered  to be responsible for the small $x$  behaviour, we show that the factorial moments  satisfy \eq{C1},
 which  reproduces the multiplicity distribution of \eq{MULTD}. This result is in agreement with the attempts \cite{GOLE}   to find solutions to the equations for the cascade of  color-singlet dipoles (see \eq{PC1}.
 
 We also suggest a procedure for finding corrections to this multiplicity distribution, which, we believe, will be useful for descriptions of the experimental data.
 
 In general, the multiplicity distribution, that has been discussed in the paper,  confirms the result of Ref.\cite{KHLE}, that the entropy  of color-singlet dipoles  is equal $S\,\,=\,\,\ln N$ in the region of small $x$, and gives the regular procedure to estimate corrections to this formula.
 
 It is worthwhile mentioning that both SCA and DA have been developed before and many technical issues that matter,  have been discussed\footnote{We thank our referee who drew our attention to Ref.\cite{BWX},  where  some technical details of our approach have been clarified.}(see, for example,  Refs.\cite{KOLEB,LRSH,LLSH,BRAUN1,BWX,LEMI} and references therein).
 
 
   
\section{Acknowledgements}
 
   We thank our colleagues at Tel Aviv university and UTFSM for
 encouraging discussions.  Our special thanks go to Michael Lublinsky  for fruitful discussion of Ref.\cite{LMM}. 
This research was supported  by 
 ANID PIA/APOYO AFB180002 (Chile) and  Fondecyt (Chile) grants  
 1180118.
 


  