High energy QCD: multiplicity dependence of quarkonia production

Abstract

In this paper we propose an approach which demonstrates the dependence of quarkonia production on the multiplicity of the accompanying hadrons. Our approach is based on the three gluon fusion mechanism, without assuming the multiplicity dependence of the saturation scale. We show, that we describe the experimental data, which has a dependence that is much steeper than the multiplicity of the hadrons.

1 Introduction

The goal of this paper is to study the multiplicity dependence of quarkonia (mainly $J/\Psi$) production in the framework of high energy QCD (see Ref. [1] for a general review). Effective QCD at high energies currently exists in two different formulations: the CGC/saturation approach  [2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18], and the BFKL Pomeron calculus  [19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46]. In this paper we restrict ourself to the BFKL Pomeron calculus, which has a more direct correspondence with the parton approach, and which provides an approximation for estimates of hadron–hadron collisions, that at present are out of the reach for the CGC approach.

Fortunately, in Refs. [44, 45] it was shown, that these two approaches are equivalent for the description of the scattering amplitude in the rapidity range :

$$\begin{aligned} Y \le \frac{2}{\Delta _{\mathrm{BFKL}}} \ln \left( \frac{1}{\Delta ^2_{\mathrm{BFKL}}}\right) \end{aligned}$$

(1)

where $\Delta _{\mathrm{BFKL}}$ denotes the intercept of the BFKL Pomeron. In this paper it is also shown, that for Eq. (1) we can use the Mueller, Patel, Salam and Iancu approximation(MPSI)  [47,48,49,50,51] for hadron–hadron scattering at high energies.

The recent experiments by ALICE [52,53,54,55,56] and STAR [57, 58], show that the cross sections for $J/\Psi$ production depends strongly on the multiplicity of accompanying hadrons. These data have stimulated theoretical discussions on the origin of such dependence (see Refs. [59,60,61,62,63]). In this paper, we develop an approach to this problem based on two ingredients. First, we assume, that the production of quarkonia stems from triple gluon fusion [61, 64, 65] (see Fig. 1). For the interaction with nuclei this mechanism is dominant  [66,67,68,69,70,71]; and it has been demonstrated in Ref. [61], that this mechanism gives a substantial contribution in hadron–hadron collisions.

Fig. 1

a The three gluon fusion mechanism of $J/\Psi$ production. b The Mueller diagram [72] for $J/\Psi$ production, which illustrates the inter-relation between three gluon fusion and the triple BFKL Pomeron interaction. The wavy lines describe the BFKL Pomerons, while the helical curves represent gluons

Second, we showed in Refs. [73, 74], 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, the multiplicity distribution has a general form:

$$\begin{aligned} \frac{\sigma _n}{\sigma _{ \mathrm in}} = \frac{1}{N} \left( \frac{N - 1}{N}\right) ^{n - 1} \end{aligned}$$

(2)

where N denotes the average number of partons.

The paper is organized as follows: in the next section we describe our approach to hadron–hadron collisions. In Sect. 3 we discuss quarkonia production in simplified Reggeon Field Theory, defined in zero transverse dimensions. In Sect. 4 we generalize the result in this toy-model approach to high energy QCD, and compare our estimates with the experimental data . We summarize our results in the Conclusions.

2 Hadron–hadron interaction in MPSI approach

This section does not contain new results, and we include it in the paper for completeness of presentation, as well as a kind of an introduction to the notation, and the main ideas.

2.1 QCD parton cascade

We start with the equation for the QCD parton cascade which can be written in the following form [1, 6, 7, 38,39,40]:

$$\begin{aligned}&\frac{\partial \,P_n\left( Y, \varvec{r }, \varvec{b};\,\varvec{r}_1,\varvec{ b}_1,\,\varvec{r}_2 , \varvec{b}_2\dots \varvec{r}_i ,\varvec{b}_i, \dots \varvec{r}_n, \varvec{b}_n \right) }{ \partial \, Y } \nonumber \\&\quad = - \sum ^n_{i=1}\,\omega _G(r_i) \, P_n\left( Y, \varvec{r }, \varvec{b};\,\varvec{r}_1,\varvec{ b}_1,\,\varvec{r}_2 , \varvec{b}_2\dots \varvec{r}_i ,\varvec{b}_i, \dots \varvec{r}_n, \varvec{b}_n \right) \nonumber \\&\qquad + \bar{\alpha }_S\,\sum ^{n-1}_{i=1} \,\frac{(\varvec{r}_i\,+\, \varvec{r}_n)^2}{(2\,\pi )\,r^2_i\,r^2_n}\, \nonumber \\&\qquad \times P_{n - 1}\left( Y, \varvec{r},\varvec{b};\,\varvec{r}_1,\varvec{b}_1, \dots (\varvec{r}_i \,+\, \varvec{r}_n), \varvec{b}_{in},\dots \varvec{r}_{n-1},\varvec{b}_n \right) \end{aligned}$$

(3)

where $P_n\left( Y ; \{r_i,b_i\}\right)$ denotes the probability to have n-dipoles of size $r_i$, at impact parameter $b_i$, and at rapidity Y.¹$\varvec{b}_{in}$ in Eq. (3) is given by $\varvec{b}_{in} \,=\,\varvec{b}_i \,+\,\frac{1}{2}\varvec{r}_i \,=\,\varvec{b}_n - \frac{1}{2}\varvec{r}_i$.

Equation (3) is a typical cascade equation in which the first term describes the depletion of the probability of n, due to one dipole decaying 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{aligned}&P_1 \left( Y = 0, \varvec{r},\varvec{b} ; \varvec{r}_1,\varvec{b}_1\right) \,\, =\,\,\,\delta ^{(2)}\left( \varvec{r} - \varvec{r}_1\right) \,\delta ^{(2)} \left( \varvec{b} - \varvec{b}_1\right) ;\nonumber \\&\quad ~~~~~~~P_{n>1}\left( Y = 0; \{r_i\}\right) \,=\,0 \end{aligned}$$

(4)

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{aligned} \sum _{n=1}^\infty \,\int \prod ^n_{i=1} d^2 r_i \,d^2 b_i \,P_n\left( Y ; \{\varvec{r}_i\,\varvec{b}_i\}\right) \,\,=\,\,1 , \end{aligned}$$

(5)

i.e. the sum of all probabilities is equal to 1.

This QCD cascade leads to the Balitsky–Kovchegov (BK) equation  [1, 8,9,10] for the amplitude, and gives the theoretical description of DIS. To see this we introduce the generating functional [6, 7]

$$\begin{aligned} Z\left( Y, \varvec{r},\varvec{b}; [u_i]\right)= & {} \sum ^{\infty }_{n=1}\int P_n\left( Y,\varvec{r},\varvec{b};\{\varvec{r}_i\,\varvec{b}_i\}\right) \nonumber \\&\times \prod ^{n}_{i=1} u\left( \varvec{r}_i\,\varvec{b}_i\right) \,d^2 r_i\,d^2 b_i \end{aligned}$$

(6)

where $u\left( \varvec{r}_i\,\varvec{b}_i\right) \equiv \,=u_i$ is an arbitrary function. The initial conditions of Eq. (4) and the sum rules of Eq. (5) require the following form for the functional Z:

$$\begin{aligned} Z\left( Y=0, \varvec{r},\varvec{b}; [u_i]\right)= & {} u\left( \varvec{r},\varvec{b}\right) ; \end{aligned}$$

(7a)

$$\begin{aligned} Z\left( Y, r,[u_i=1]\right)= & {} 1; \end{aligned}$$

(7b)

Multiplying both terms of Eq. (3) by $\prod ^{n}_{i=1} u\left( \varvec{r}_i \,\varvec{b}_i\right)$ and integrating over $r_i$ and $b_i$, we obtain the following linear functional equation [39, 40];

$$\begin{aligned}&\frac{\partial Z\left( Y, \varvec{r},\varvec{b}; [u_i]\right) }{\partial \,Y} \nonumber \\&\quad =\int d^2 r'\, K\left( \varvec{r}',\varvec{r} - \varvec{r'}|\varvec{r}\right) \Bigg ( - u\left( r, b\right) \nonumber \\&\qquad + u\left( \varvec{r}',\varvec{b} + \frac{1}{2}(\varvec{r} - \varvec{r}') \right) \,u\left( \varvec{r} - \varvec{r}',\varvec{b} - \frac{1}{2}\varvec{r}'\right) \Bigg ) \frac{\delta \,Z}{\delta \,u\left( r, b \right) }; \end{aligned}$$

(8a)

$$\begin{aligned}&K\left( \varvec{r}',\varvec{r} - \varvec{r'}|\varvec{r}\right) \,=\frac{\bar{\alpha }_S}{2 \,\pi }\frac{r^2}{r'^2\,(\varvec{r} - \varvec{r}')^2} ;\nonumber \\&\omega _G\left( r\right) \,\,=\,\,\int d^2 r' K\left( \varvec{r}',\varvec{r} - \varvec{r'}|\varvec{r}\right) ; \end{aligned}$$

(8b)

Searching for a solution of the form $Z\left( [ u(r_i,b_i,Y)]\right)$ for the initial conditions of Eq. (7a), Eq. (8a) can be re-written as the non-linear equation  [6, 7]:

$$\begin{aligned} \frac{\partial Z\left( Y, \varvec{r},\varvec{b}; [u_i]\right) }{\partial \,Y}= & {} \int d^2 r' K\left( \varvec{r}',\varvec{r} - \varvec{r'}|\varvec{r}\right) \nonumber \\&\times \Bigg \{Z\left( r', \varvec{b} + \frac{1}{2}(\varvec{r} - \varvec{r}'); [u_i]\right) \nonumber \\&\times Z\left( \varvec{r} - \varvec{r'}, \varvec{b} - \frac{1}{2}\varvec{r}'; [u_i]\right) \, - \,Z\left( Y, \varvec{r},\varvec{b}; [u_i]\right) \Bigg \}\nonumber \\ \end{aligned}$$

(9)

Therefore, the QCD parton cascade of Eq. (3) takes into account non-linear evolution. Generally speaking the scattering amplitude can be written in the form [10, 39, 40]:

$$\begin{aligned} N(Y,\,r,\,b)= & {} - \,\sum ^{\infty }_{n=1}\,(-1)^n\, \rho ^p_n(r_1,\, b_1,\,\ldots \,r_n,\,b_n\,;Y - Y_0) \nonumber \\&\times \prod ^n_{i =1}\,N(Y_0,\,r_i,\,b_i)\,\, d^2\, r_i\, \, d^2\, b_i \,. \end{aligned}$$

(10)

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 $\rho ^p_n(r_1, b_1,\ldots \,,r_n, b_n)$ are defined as follows:

$$\begin{aligned}&\rho ^p_n(r_1, b_1\, \ldots \,,r_n, b_n; Y - Y_0)\,=\,\frac{1}{n!}\,\nonumber \\&\quad \times \prod ^n_{i =1} \,\frac{\delta }{\delta u_i } \,Z\left( Y - Y_0;\,[u] \right) |_{u=1} \end{aligned}$$

(11)

For $\rho _n$ we obtain [39, 40] :

$$\begin{aligned}&\frac{\partial \,\rho ^p_n(r_1, b_1\,\ldots \,,r_n, b_n)}{ \bar{\alpha }_s\,\partial \,Y}\nonumber \\&\quad = - \sum _{i=1}^n \,\,\omega (r_i)\,\,\rho ^p_n(r_1, b_1\,\ldots \,,r_n, b_n) \nonumber \\&\qquad + 2\,\sum _{i=1}^n\, \int \,\frac{d^2\,r'}{2\,\pi }\, \frac{r'^2}{r^2_i\,(\varvec{r}_i - \varvec{r}')^2}\, \rho ^p_n(\ldots \,r', b_i-r'/2\dots )\nonumber \\&\qquad \, +\,\sum _{i=1}^{n-1}\,\frac{(\varvec{r}_i + \varvec{r}_n)^2}{(2\,\pi )\,r^2_i\,r^2_n}\, \rho ^p_{n-1}(\ldots \,(\varvec{r}_i\,+\,\varvec{r}_n), b_{in}\dots ).\nonumber \\ \end{aligned}$$

(12)

For $\rho _1$ we have the linear BFKL equation [19,20,21]:

$$\begin{aligned} \frac{\partial \,\rho ^p_1(Y; r_1, b)}{ \bar{\alpha }_S\,\partial \,Y}= & {} - \,\omega _G\left( r_1\right) \rho ^p_1(Y; r_1, b) \nonumber \\&+ 2\,\int \,\frac{d^2\,r'}{2\,\pi }\, \frac{r'^2}{r^2_1\,(\varvec{r}_1 - \varvec{r}')^2}\, \bar{\rho }^p_1\left( Y, r',b\right) \nonumber \\ \end{aligned}$$

(13)

However, to obtain the BK equation for the scattering amplitude we need to use Eq. (10), in which we introduce the amplitude of interaction of the dipole with the target at low energies. Using Eqs. (8a),  (10) and (11) , we can obtain the non-linear BK equation from Eq. (9) in the following form [10]

$$\begin{aligned}&\frac{\partial }{\partial Y} N\left( \varvec{r}, \varvec{b} , Y \right) \nonumber \\&\quad =\int d^2 r'\,K\left( \varvec{r}', \varvec{r} - \varvec{r}'| \varvec{r}\right) \Bigg \{N\left( \varvec{r}',\varvec{b} - \frac{1}{2}\left( \varvec{r} - \varvec{r}' \right) , Y\right) \nonumber \\&\qquad + N\left( \varvec{r} - \varvec{r}', \varvec{b} - \frac{1}{2}\varvec{r}', Y\right) \,\,- \,\,N\left( \varvec{r},\varvec{b},Y \right) \nonumber \\&\qquad {-} N\left( \varvec{r} {-} \varvec{r}', \varvec{b} - \frac{1}{2}\varvec{r}', Y\right) \, N\left( \varvec{r}',\varvec{b} {-} \frac{1}{2}\left( \varvec{r} {-} \varvec{r}' \right) , Y\right) \Bigg \}\nonumber \\ \end{aligned}$$

(14)

2.2 The interaction of two dipoles at high energies

We first consider the simplest case of scattering, the high energy interactions of two dipoles with sizes r and R and with $r \,\sim \, R$. In Refs. [44, 45] it is shown that in the limited range of rapidities, given by Eq. (1), we can safely apply the Mueller, Patel, Salam and Iancu approach for this scattering  [47,48,49,50,51] (see Fig. 2a).

Fig. 2

Scattering amplitude for the interaction of two dipoles with sizes: r and R at high energy in the MPSI approach (see a and b). The amplitudes of interaction of two dipoles in the Born approximation of perturbative QCD ( $N\left( r_i,r'_i,b''_i\right)$ in Eq. (15)) are shown as white circles. The wavy lines denote the BFKL Pomerons. c Shows the Mueller diagram [72] for inclusive production of gluons

The scattering amplitude in this approach can be written in the following form [39, 40]:

$$\begin{aligned}&N\left( Y, r, R, b \right) \nonumber \\&\quad =\, - \, \sum ^{\infty }_{n =1}\,\,(-1)^n\, \int \,\,\rho ^t_n\left( \varvec{r}_1, \,\varvec{b}'_1 ,\ldots \,, \varvec{r}_n,\, \varvec{b}'_n;\,\frac{1}{2}Y\right) \,\,\nonumber \\&\qquad \times \rho ^p_n\left( \varvec{r}'_1, \, \varvec{b} - \varvec{b}'_1 - \varvec{b}_1''\,,\ldots \,,\varvec{r}_n,\, \varvec{b} - \varvec{b}'_n - \varvec{b}_n'';\, -\frac{1}{2}Y\right) \nonumber \\&\qquad \times \,\prod ^n_{i =1}\, \, d^2\, r_i \, \,\,\prod ^n_{j =1}\, \,d^2 \,r'_j \, d^2 b'_j d^2 b''_j\,N^\mathrm{BA}\left( r_i, r'_i,b''_{i}\right) \end{aligned}$$

(15)

where $\rho ^t_n$ and $\rho ^p_n$ denote the parton densities in the target and projectile, respectively. These densities can be calculated from $P_n$ using Eq. (11). $N^\mathrm{BA}$ is the scattering amplitude of two dipoles in the Born approximation of perturbative QCD (see Fig. 2). Equation (15) simply states that we can consider the QCD parton cascade of Eq. (3) generated by the dipole of size r for the c.m.f. rapidities from 0 to $\frac{1}{2}Y$, and the same cascade for the dipole of the size R, for the rapidities from 0 to $-\frac{1}{2}Y$. One can see that Eq. (15) is the t-channel unitarity re-written in a form, convenient for applying the evolution of the parton cascade in the form of Eq. (12).

Generally speaking, for a dense system of partons at $Y \, =\,0$ n-dipoles from upper cascade could interact with m dipoles from the lower cascade, with the amplitude $N^m_n\left( \{r_i\},\{r'_j\}\right)$ [39, 40]. In Eq. (15) we assume that the system of dipoles that has been created at $Y =0$ is not very dense, at least for the range of rapidities given by Eq. (1). In this case

$$\begin{aligned} N^m_n\left( \{r_i\},\{r'_j\}\right) \,\,=\,\,\delta _{n,m}\prod ^n_{j =1}\,\left( - 1\right) ^{n - 1}N^\mathrm{BA}\left( r_i, r'_i,b''_i\right) \nonumber \\ \end{aligned}$$

(16)

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. 2b.

The average number of dipoles at $Y = 0$ is determined by the inclusive cross section, which is given by the diagram of Fig. 2c and which can be written at $y \rightarrow 0$ as follows [75]:

$$\begin{aligned} \frac{d \sigma }{d y \,d^2 p_{T}}= & {} \frac{2C_F}{\alpha _s (2\pi )^4}\,\frac{1}{p^2_T} \int d^2 \varvec{r}_T\,e^{i \varvec{p}_T\cdot \varvec{r}_T}\!\nonumber \\&\times \int \!\! d^2 b\,\nabla ^2_T\,N^\mathrm{\mathrm{BFKL}}\left( \frac{1}{2}Y ;r, r_T; b \right) \nonumber \\&\times \int \! d^2 B \,\nabla ^2_T\,N^\mathrm{\mathrm{BFKL}}\left( y_2 = -\frac{1}{2}Y ;R, r_T; B \right) \nonumber \\ \end{aligned}$$

(17)

The average number of dipoles that enters the multiplicity distribution of Eq. (27), is equal $\bar{n} \,=\,N\,=\int \frac{d^2 p_T}{(2 \pi )^2} \frac{d \sigma }{d y \,d^2 p_{T}}\Big / \sigma _{in} \,\propto \,\exp \left( \Delta _\mathrm{BFKL}\,Y\right)$² only if we assume that $\sigma _{in} \,\sim \mathrm{Const}$. Indeed, the enhanced diagrams of Fig. 2b lead to the inelastic cross section which is constant at high energy.

2.3 Hadron–hadron collisions

In this paper we view a hadron as a dilute system of dipoles and use Eq. (17) for the average multiplicity, together with the multiplicity distribution of Eq. (1). In particular, we assume that Eq. (16) is correct, and the system of partons that is produced at c.m. rapidity $y^*$=0 is a dilute system. However, we are aware that Eq. (17) does not describe the experimental increase of the average multiplicity, which from Eq. (17) is $\bar{n} \,\propto \,\exp \left( \Delta _\mathrm{BFKL}\,Y\right)$. The experimental data can be described in the framework of the CGC/saturation approach in which $N^\mathrm{\mathrm{BFKL}}$ were replaced by $N^{\mathrm{BK}}$ [76]. Hence, we cannot view hadrons as a dilute system of dipoles, but rather have to consider them as a dense system of dipoles. For such a situation we expect that $\bar{n} \,\propto \,Q^2_s(Y)/\bar{\alpha }_S$ (see Refs. [1, 76,77,78,79,80,81,82].

3 Reggeon field theory in zero transverse dimensions

3.1 Multiplicity distribution: a recap

In the parton model [83,84,85,86] 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 a dipole of fixed size. In this case the cascade equation (see Eq. (3)) takes the following simple form:

$$\begin{aligned} \frac{d P_n\left( Y\right) }{d Y}\,\,=\,\,- \Delta \,n\, P_n\left( Y\right) \,\,+\,\,\left( n - 1 \right) \Delta P_{n-1}\left( Y .\right) \end{aligned}$$

(18)

where $P_n\left( Y \right)$ denotes the probability to find n dipoles (of a fixed size in our model) at rapidity Y, and $\Delta$ denotes the intercept of the BFKL Pomeron.

Instead of the generating functional of Eq. (6), we can introduce the generating function:

$$\begin{aligned} Z\left( Y, u\right) \,\,=\,\,\sum _n \,P_n\left( Y\right) \,u^n \end{aligned}$$

(19)

where u are numbers.

At the initial rapidity $Y=0$, we have only one dipole, so $P_1\left( Y = 0 \right) =1$ and $P_{n > 1} \,=\,0$ (so the state is only one dipole); at $u = 1$, $Z\left( Y, u=1\right) \,\,=\,\,\sum _n P\left( y \right) \,=\,1$. These two properties determine the initial and the boundary conditions for the generating function which simplify Eqs. (7a) and  (7b)

$$\begin{aligned} Z\left( Y = 0,u\right) \,\,=\,\,u;~~~~~~~~~Z\left( Y, u=1\right) \,\,=\,\,1. \end{aligned}$$

(20)

Equation (18) takes the following form for the generating function:

$$\begin{aligned} \frac{\partial Z\left( Y, u\right) }{ \partial Y}\,\,=\,\,- \Delta \, u \left( 1 - u\right) \frac{\partial Z\left( Y, u\right) }{\partial u} . \end{aligned}$$

(21)

The general solution to Eq. (21) is an arbitrary function ($Z\left( z \right)$) of the new variable: $z\,\,=\,\,\Delta \,Y\,\,+\,\,f(u)$, with f(u) from the following equation:

$$\begin{aligned} 1\,\,=\, - u\,\left( 1 - u\right) \,f'_u\left( u\right) ~~~~~f\left( u\right) \,\,=\,\,\ln \left( \frac{u - 1}{u}\right) \,\,+\,\,C_1 \end{aligned}$$

(22)

The form of arbitrary function stems from the initial condition of Eq. (20)

$$\begin{aligned} Z\left( z\left( Y=0\right) \right) \,=\,u; \end{aligned}$$

(23)

Since $u\,\,=\,\,1\Big /\left( 1 - e^z\right)$ we obtain that

$$\begin{aligned} Z\left( Y,\,u\right) \,\,= & {} \,\,\frac{u\,\,e^{ - \Delta \, Y}}{1\,\,+\,\,u\,\, \left( e^{ - \Delta \,Y} - 1\right) }\,\,\nonumber \\= & {} \,\,u\,e^{ - \Delta \, Y}\,\sum ^\infty _{n=1}\,u^n \left( 1 - e^{ - \Delta \,Y}\right) ^{n - 1} . \end{aligned}$$

(24)

Note, that $Z\left( Y,u=1\right) = 1$, as it should be from Eq. (20).

On the other hand, we can re-write Eq. (21) in the form of the non-linear equation using Eq. (22): viz.

$$\begin{aligned} \,\,\frac{\partial \,Z}{\partial \,Y}\,\,=\, - \Delta \left( Z \, - \, Z^2\right) . \end{aligned}$$

(25)

Comparing Eq. (24) with Eq. (19) one can see that

$$\begin{aligned} P_n\left( Y \right) \,\,=\,\,e^{ - \Delta \, Y }\left( 1 - e^{ - \Delta \,Y}\right) ^{n - 1} . \end{aligned}$$

(26)

Since from Eq. (26) it follows that the average $n \,=\,N$ is equal to $N\,= \exp \left( \Delta \,Y\right)$ Eq. (26) can be re-written in the form of Eq. (2):

$$\begin{aligned} P_n\left( N\right) \,\,=\,\,\frac{1}{N}\,\Big ( 1 - \frac{1}{N}\Big )^{n - 1} \end{aligned}$$

(27)

Fig. 3

Mueller diagrams [72] for inclusive $J/\Psi$ production in the hadron–hadron collisions. The wavy lines denote the Pomeron Green’s functions

3.2 Quarkonia production

As we have discussed in the introduction, we assume the production of heavy quakonia stems from three gluon fusion (see Fig. 1), and it is intimately related to the triple Pomeron interaction. The Mueller diagrams for inclusive $J/\Psi$ production are shown in Fig. 3, where the wavy lines denote the Pomeron Green’s function ($G_{I\!\!P}$) which is equal to

$$\begin{aligned} G_{I\!\!P}\left( Y\right) \,\,=\,\,e^{\Delta \,Y}. \end{aligned}$$

(28)

The MPSI approach for inclusive production with fixed multiplicity of produced hadrons is shown in Fig. 4.

Therefore, from this figure we see that the structure of the parton cascade for the quarkonia production is quite different. In particular, for the part of the events whose weight is determined by the contribution of Fig. 3b, the initial conditions for the parton cascade is not the ones of Eq. (20) however, they have the form:

$$\begin{aligned} Z\left( Y = 0,u\right) \,\,=\,\,u^2;~~~~~~~~~Z\left( Y, u=1\right) \,\,=\,\,1. \end{aligned}$$

(29)

This means, that for this cascade we need to find the arbitrary function $Z\left( z\right)$ from the following equation

$$\begin{aligned} Z\left( z\left( Y=0\right) \right) \,=\,u^2; \end{aligned}$$

(30)

The solution is

$$\begin{aligned} Z\left( Y, u\right)= & {} \frac{u^2\,e^{- 2\,\Delta \,Y}}{\left( 1 \,\,+\,\,u\left( e^{ - \,\Delta \,Y}\, - \,1\right) \right) ^2}\nonumber \\= & {} u^2\,e^{ - 2\,\Delta \, Y}\,\sum ^\infty _{n=1}\left( n - 1\right) \,u^n \left( 1 - e^{ - \Delta \,Y}\right) ^{n -2} \end{aligned}$$

(31)

Fig. 4

The example of Mueller diagrams [72] for inclusive $J/\Psi$ production in hadron–hadron collisions with fixed multiplicity of produced hadrons and the parton cascade, which they describe. The exchange of Pomerons do not cancel each other due to AGK cutting rules [87], since we fixed the multiplicity in the final state. The wavy lines denote the Pomeron Green’s functions. The black circles indicate the triple Pomeron vertices. The solid lines correspond to partons. $\Delta y$ is the rapidity window in which the multiplicity of the soft hadron is measured. This window is situated in central rapidity region with rapidity about $\frac{1}{2}Y$

Equation (31) leads to a different multiplicity distribution in comparison with Eq. (27): viz.

$$\begin{aligned} P^{(2)}_n\,\,=\,\,\frac{1}{N^2} \left( n\, - \,1\right) \left( 1 - \frac{1}{N}\right) ^{n - 2} \end{aligned}$$

(32)

Finally, the cross section for quarkonia production with given multiplicity (n) of produced hadrons in the rapidity window $\Delta y$ (see Fig. 4), is equal to

$$\begin{aligned} \frac{d \sigma _n^{\mathrm{J/\Psi }}}{ d y}= & {} \frac{d \sigma ^\mathrm{J/\Psi }_\mathrm{incl}}{ d y}\left( Fig.~3a\right) \frac{n}{\langle n^{(1)}\rangle } \, \, P^{(1)}_n\left( N\right) \nonumber \\&+ \frac{d \sigma ^\mathrm{J/\Psi }_\mathrm{incl}}{ d y}\left( Fig.~3b\right) \frac{n}{\langle n^{(2)}\rangle } \, \, P^{(2)}_n\left( N\right) \end{aligned}$$

(33)

where $\langle n^{(1)}\rangle$ and $\langle n^{(2)}\rangle$ are average multiplicities of the produced hadrons in the parton cascades which are shown in Fig. 4a and in Fig. 4b, respectively.

In Eq. (33) the first and the second terms correspond to parton cascades that are shown in Fig. 4a and in Fig. 4b, respectively. The appearance of the factor $n/\langle n^{(i)}\rangle$ , which is the number of parton ladders, stems from the fact that $J/\Psi$ can be produced from every parton ladder (cut Pomeron) [60, 62, 63]. It should be stressed that the number of ladders at rapidity y from which the $J/\Psi$ is produced, is the same as the number of ladders from which the soft hadrons are produced in the MPSI approach [47,48,49,50,51]. It follows from the fact that integration over rapidities ($y'$) of the triple Pomeron vertices in Fig. 4 leads to $Y - y' \propto 1/\Delta$ for Pomerons in upper part of the diagram, and to $y' \propto 1/\Delta$ in the lower part of the diagram. $\Delta$ denotes the Pomeron intercept.

Summing Eq. (33) we obtain that

$$\begin{aligned} \sum _n \frac{d \sigma _n^\mathrm{J/\Psi }}{ d y}= & {} \frac{d \sigma ^\mathrm{J/\Psi }_\mathrm{incl}}{ d y}\left( Fig.~3a\right) \nonumber \\&+\,\, \frac{d \sigma ^\mathrm{J/\Psi }_\mathrm{incl}}{ d y}\left( Fig.~3b\right) \,\,=\,\,\frac{d \sigma ^\mathrm{J/\Psi }_\mathrm{incl}}{ d y} \end{aligned}$$

(34)

which coincide with the Mueller diagram approach, shown in Fig. 3. Introducing $\kappa \,\, = \frac{d \sigma ^\mathrm{J/\Psi }_\mathrm{incl}}{ d y}\left( Fig.~3b\right) \Bigg / \frac{d \sigma ^\mathrm{J/\Psi }_\mathrm{incl}}{ d y}\left( Fig.~3a\right)$ we re-write Eq. (33) in the following form:

$$\begin{aligned} \frac{\frac{d \sigma _n^\mathrm{J/\Psi }}{ d y}}{\frac{d \sigma ^\mathrm{J/\Psi }_\mathrm{incl}}{ d y}}= & {} \frac{1}{1\,\,+\,\,\kappa }\Bigg ( \frac{n}{\langle n^{(1)}\rangle } \, \, P^{(1)}_n\left( N\right) \nonumber \\&+\,\,\, \kappa \frac{n}{\langle n^{(2)}\rangle } \, \, P^{(2)}_n\left( N\right) \Bigg ) \end{aligned}$$

(35)

where $\kappa$ is equal to [61]

$$\begin{aligned} \kappa \,\,=\,\,e^{ 2\,\Delta \left( \frac{1}{2}Y \, - \,y\right) } \end{aligned}$$

(36)

In Eq. (33) $P^{(1)}_n\left( N\right) \,\,\equiv \,\,P_n\left( N\right)$ of Eq. (27).

The cross section of produced gluons (hadrons) is proportional to

$$\begin{aligned} \frac{d \sigma ^\mathrm{prod. gl.}_n}{ d y} \,\,=\,\, \, \, \frac{d \sigma ^\mathrm{prod. gl.}_\mathrm{incl}}{ d y} P^{(1)}_n\left( N\right) \end{aligned}$$

(37)

From Fig. 5, one can see that $P^{(2)}_n\left( N\right)$ and $P^{(1)}_n\left( N\right)$ have different dependance on $z = n/N$. The average number of gluons is chosen to be the mean multiplicity of hadrons in the rapidity window $|\eta |\,\le \, 0.9$ measured at W=13TeV.

Fig. 5

a $P^{(2)}_n\left( N\right)$ and $P^{(1)}_n\left( N\right)$ versus $z = \frac{n}{N}$. N is taken to be equal 10. b The ratio of $P^{(2)}_n\left( N\right) /P^{(1)}_n\left( N\right)$ for $N=10$ and for large N (see Eq. (39))

Therefore, $\frac{d \sigma ^\mathrm{J/\Psi }}{ d y}$ is not proportional to $\frac{d \sigma ^\mathrm{prod.\, gl.i}}{ d y}$, but shows non-linear dependance, which we will discuss below. It should be stressed that such dependance stems from triple Pomeron mechanism of quarkonia production, as noted in Refs. [61, 63]. At large N we have

$$\begin{aligned}&P^{(2)}_n\,\,\xrightarrow {N\,\gg \,1}\,\,\frac{1}{N} z\,e^{-z};\nonumber \\&P^{(1)}_n\,\,\xrightarrow {N\,\gg \,1}\,\,\frac{1}{N} \,e^{-z}; \end{aligned}$$

(38)

leading to

$$\begin{aligned} \frac{ P^{(2)}_n}{P^{(1)}_n} \,\,\xrightarrow {N\,\gg \,1}\,\,z; \end{aligned}$$

(39)

It is worth mentioning that the average multiplicity of $P^{(2)}$ distribution is equal to

$$\begin{aligned} < n^{(2)} >\,\,=\,\,\sum ^\infty _{n=2}\,n\, P^{(2)}_n\left( N\right) \,\,=\,\,2\,N \end{aligned}$$

(40)

Hence, for the multiplicity distribution $P^{(2)}$ the average number of accompanying gluons(hadrons) is twice larger than in the distribution $P^{(1)}$ . Equation (40) means that the ratio $\frac{n^{(2)}}{ < n^{(2)} > }\,\,=\,\,\frac{n}{ 2\,N}$.

4 Structure of QCD parton cascade

4.1 Multiplicity distribution

In QCD, to find the multiplicity distribution for hadron–hadron scattering in QCD using MPSI approach [47,48,49,50,51], we need to evaluate (see Eq. (6))

$$\begin{aligned} \tilde{P}_n\left( Y, r\right) \,\,=\,\, \int P_n\left( Y,\varvec{r},\varvec{b};\{\varvec{r}_i\,\varvec{b}_i\}\right) \prod ^{n}_{i=1} \,d^2 r_i\,d^2 b_i \end{aligned}$$

(41)

However, in the MPSI approach it is more natural to introduce moments (see Eq. (12)):

$$\begin{aligned} M^p_n\left( Y, r\right)= & {} \int _{r} \prod ^n_{i=1} d^2 r_i d^2 b \,\rho ^p_n\left( Y, \{ r_i\},b\right) \nonumber \\= & {} \int _{r} \prod ^n_{i=1}\frac{ d^2 r_i }{r^2_i } d^2 b\, \bar{\rho }^p_n\left( Y, \{r_i\},b\right) \end{aligned}$$

(42)

The integration over $r_i$ depends on the size of the initial dipole, which generates the cascade. In DIS the natural integration stems from $r_i \,>\,r$.

4.1.1 Several first iterations.

We start from the first several iteration of Eq. (12), which can be re-written for $\,\bar{\rho }^p_n(r_1, b_1\,\ldots \,,r_n, b_n)$ in the form

$$\begin{aligned}&\frac{\partial \,\bar{\rho }^p_n(r_1, b_1\,\ldots \,,r_n, b_n)}{ \bar{\alpha }_s\,\partial \,Y} \nonumber \\&\quad = -\,\sum _{i=1}^n \,\,\omega (r_i)\,\,\bar{\rho }^p_n(r_1, b_1\,\ldots \,,r_n, b_n)\nonumber \\&\qquad +2\,\sum _{i=1}^n\, \int \,\frac{d^2\,r'}{2\,\pi }\, \frac{1}{(\varvec{r}_i - \varvec{r}')^2}\, \bar{\rho }^p_n(\ldots \,r', b_i-r'/2\dots )\nonumber \\&\qquad +\,\sum _{i=1}^{n-1}\, \bar{\rho }^p_{n-1}(\ldots \,(\varvec{r}_i\,+\,\varvec{r}_n), b_{in}\dots ). \end{aligned}$$

(43)

For the first iteration $\bar{\rho }_1\left( Y; r_1,b_1\right)$, we obtain the BFKL equation:

$$\begin{aligned} \frac{\partial \,\bar{\rho }^p_1(Y; r_1, b)}{ \bar{\alpha }_S\,\partial \,Y}= & {} - \,\omega _G\left( r_1\right) \bar{\rho }^p_1(Y; r_1, b)\nonumber \\&+\,2\,\int \,\frac{d^2\,r'}{2\,\pi }\, \frac{1}{(\varvec{r}_1 - \varvec{r}')^2}\, \bar{\rho }^p_1\left( Y, r',b\right) \nonumber \\= & {} \int d^2 \,r' K\left( r_1,r'\right) \bar{\rho }^p_1\left( Y, r',b\right) \end{aligned}$$

(44)

with the solution

$$\begin{aligned}&\bar{\rho }^p_1(Y; r_1, b)\,\,=\,\,\int ^{\epsilon \,+\,i\,\infty }_{\epsilon - i\,\infty } \frac{d \omega }{ 2\,\pi \,i} \int ^{\epsilon \,+\,i\,\infty }_{\epsilon - i\,\infty } \nonumber \\&\quad \frac{d \gamma }{ 2\,\pi \,i} \,e^{ \omega \,Y\,\,+\,\,\gamma \,\xi _1}\frac{1}{\omega - \bar{\alpha }_S\chi \left( \gamma \right) } \tilde{\rho }^p_{in,1}(\gamma , b) \end{aligned}$$

(45)

where $\xi _1\,\,=\,\,\ln \left( r^2_1\Lambda _\mathrm{QCD}^2\right)$ and

$$\begin{aligned}&\omega \left( \gamma \right) =\bar{\alpha }_S\,\chi \left( \gamma \right) =\bar{\alpha }_S\left( 2 \psi \left( 1\right) \right. \left. - \psi \left( \gamma \right) - \psi \left( 1 - \gamma \right) \right) \nonumber \\&\hbox {d.a} \,\xrightarrow {\gamma \rightarrow \frac{1}{2}} \,\omega _0\,\,+\,\,D\,\left( \gamma - \frac{1}{2}\right) ^2 +\,\,\mathcal{O}\left( \left( \gamma - \frac{1}{2}\right) ^3\right) \,\,\nonumber \\= & {} \,\,\bar{\alpha }_S4 \ln 2 \,\,+\,\,\bar{\alpha }_S14 \zeta \left( 3\right) \left( \gamma - \frac{1}{2}\right) ^2 \,\,+\,\,\mathcal{O}\left( \left( \gamma - \frac{1}{2}\right) ^3\right) \nonumber \\ \end{aligned}$$

(46)

where d.a denotes ‘diffusion approximation’ and $\psi (z)$ is Euler gamma function (see [88] formula 8.36). $\tilde{\rho }^p_{in,1}(\gamma , b)$ has to be found from the initial conditions. From Eq. (45) we can obtain $M^p_1\left( Y,r,b\right)$ (see Eq. (42)) which has the following form:

$$\begin{aligned} M^p_1\left( Y,r\right)= & {} \int ^{\epsilon \,+\,i\,\infty }_{\epsilon - i\,\infty } \frac{d \omega }{ 2\,\pi \,i} \int ^{\epsilon \,+\,i\,\infty }_{\epsilon - i\,\infty } \frac{d \gamma }{ 2\,\pi \,i} \nonumber \\&\times e^{ \omega \,Y\,\,+\,\,\gamma \,\xi }\frac{1}{\omega - \bar{\alpha }_S\chi \left( \gamma \right) } M^p_{in,1}(\gamma ) \end{aligned}$$

(47)

which satisfies the following equation:

$$\begin{aligned}&\frac{ \partial \, M^p_1\left( Y, r\right) }{\partial \,\bar{\alpha }_S\,Y}\,\,=\,\,\,\,\,\int d^2 \,r' K\left( r,r'\right) M^p_1\left( Y, r'\right) \,\,\nonumber \\&\quad \xrightarrow {r'\,\gg \,r}\,\,\,\int _r\frac{d \,r'^2}{r'^2} \,\,M^p_1\left( Y, r'\right) \end{aligned}$$

(48)

The equation for the next iteration: $\bar{\rho }_2$, takes the form: equation for $\rho ^p_2$ can be re-written in the following form for $\bar{\rho }^p_2$:

$$\begin{aligned}&\frac{\partial \,\bar{\rho }^p_2(Y; r_1, r_2, b)}{ \bar{\alpha }_S\,\partial \,Y}\nonumber \\&\quad =\int d^2\,r'\,K\left( r_1, r'\right) \bar{\rho }^p_2\left( Y, r',b, r_2,b\right) \nonumber \\&\qquad +\,\,\,\int d^2\,r'\,K\left( r_2, r'\right) \bar{\rho }^p_2\left( Y, r_1, r', b\right) \nonumber \\&\qquad +\,\,\bar{\rho }^p_{1}\left( Y; \varvec{r}_1\,+\,\varvec{r}_2, b\right) \end{aligned}$$

(49a)

$$\begin{aligned}&\quad \xrightarrow {r'\,\gg \,r_i} \int _{r_1} \frac{d r'^2}{r'^2} \bar{\rho }^p_2\left( Y, r',b, r_2,b\right) \nonumber \\&\qquad + \int _{r_2} \frac{d r'^2}{r'^2} \bar{\rho }^p_2\left( Y, r_1,b, r',b\right) \,\,+\,\,\bar{\rho }^p_{1}\left( Y; \varvec{r}_1\,+\,\varvec{r}_2, b\right) \nonumber \\ \end{aligned}$$

(49b)

For simplicity we re-write Eq. (49a) in the log approximation following Ref. [73] (see Eq. (49b)). Rewriting Eq. (49b) for $M^p_2$ we obtain:

$$\begin{aligned}&\frac{\partial \,M^p_2(Y; r)}{ \bar{\alpha }_S\,\partial \,Y}\nonumber \\&\quad =\,\, \int _r \frac{d r'^2}{r'^2} \,\,\Bigg \{ 2\,M^p_2\left( Y, r' \right) \,\,+\,\,M^p_1\left( Y, r'\right) \Bigg \}\nonumber \\&\quad =\,\, 2\,\int _r \frac{d r'^2}{r'^2} \,\,\,M^p_2\left( Y, r' \right) \,\,+\,\,\frac{\partial }{ \partial \,\bar{\alpha }_S\,Y} M^p_1\left( Y, r \right) \end{aligned}$$

(50)

In the last term of Eq. (50) we used Eq. (48). The solution to Eq. (50) takes the form:

$$\begin{aligned} M^p_2\left( Y, r\right)= & {} \int ^{\epsilon \,+\,i\,\infty }_{\epsilon - i\,\infty } \frac{d \omega }{ 2\,\pi \,i} \int ^{\epsilon \,+\,i\,\infty }_{\epsilon - i\,\infty } \frac{d \gamma }{ 2\,\pi \,i} \,e^{ \omega \,Y\,\,+\,\,\gamma \,\xi }\nonumber \\&\times \frac{\bar{\alpha }_S\chi \left( \gamma \right) }{\left( \omega \, - \,2\,\bar{\alpha }_S\chi \left( \gamma \right) \right) \,\left( \omega \, - \,\bar{\alpha }_S\chi \left( \gamma \right) \right) } \end{aligned}$$

(51)

with $\chi \left( \gamma \right) = 1/\gamma$. One can check that $M^p_2\left( Y, r\right) \,\xrightarrow {Y \,\rightarrow \,0}\,\,0$, which is the correct initial condition for one dipole of size r at $Y=0$ , which generates the parton cascade. Actually, Eq. (50) describes $M^p_2\left( Y, r\right)$ for the full BFKL kernel. Indeed, considering Eq. (49a) we can integrate this equation over $r_1$ and $r_2$, to obtain the equation for $M^p_2$. The last term has the following form

$$\begin{aligned}&\int _r \frac{d^2 r_1}{2\,\pi } \frac{1}{r^2_1} \frac{d^2 r_2}{2\,\pi }\int d^2 b \frac{1}{r^2_2} \bar{\rho }(Y,r_{12}, b)\nonumber \\&\quad = \int _r \frac{d^2 r_1}{2\,\pi } \frac{1}{r^2_1} \frac{d^2 r_2}{2\,\pi }\int d^2 b \, \frac{1}{\left( \varvec{r}_1 - \varvec{r}_{12}\right) ^2} \,\bar{\rho }(Y,r_{12}, b) \end{aligned}$$

(52)

Note that $r^2_2\,=\,\left( \varvec{r}_1 - \varvec{r}_{12}\right) ^2$ can never approach zero, since $r_2> r$. Removing this restriction, we can re-write

$$\begin{aligned}&\int _r \frac{d^2 r_2}{2\,\pi }\, \frac{1}{\left( \varvec{r}_1 - \varvec{r}_{12}\right) ^2} \,\bar{\rho }(Y,r_{12}, b)\nonumber \\&\quad = \int _0 \frac{d^2 r_2}{2\,\pi }\, \frac{1}{\left( \varvec{r}_1 - \varvec{r}_{12}\right) ^2} \,\bar{\rho }(Y,r_{12}, b)\, - \,\ln r^2_1\, \bar{\rho }(Y,r_{1}, b)\nonumber \\ \end{aligned}$$

(53)

The reggeization term in Eq. (53) describes the contribution of $r_2\,\rightarrow \,0$. Plugging Eq. (53) in the last term of Eq. (49a) integrated over $r_1$ and $r_2$, one can see that we reproduce Eq. (50) for the full BFKL kernel. Hence Eq. (51) is the solution with $\chi \left( \gamma \right)$ which is given by Eq. (46).

4.1.2 General solution

The equation for $M^p_n\left( Y, r\right)$ has the following general form:

$$\begin{aligned}&\frac{\partial \,M^p_n(Y; r)}{ \bar{\alpha }_S\,\partial \,Y}\,\,\nonumber \\&\quad =\,\, \int d^2\,r'\,K\left( r, r'\right) \,\,\Bigg \{ n\,M^p_2\left( Y, r' \right) \,\,+\,\,(n - 1)M^p_1\left( Y, r'\right) \Bigg \}\nonumber \\&\quad =\,\, n\,\int d^2\,r'\,K\left( r, r'\right) \,\,\,M^p_2\left( Y, r' \right) \,\,+\,\,\left( n - 1\right) \frac{\partial }{ \partial \,\bar{\alpha }_S\,Y} M^p_1\left( Y, r \right) \nonumber \\ \end{aligned}$$

(54)

The solution to this equation [73], which gives $M^p_1\left( Y=0, r\right) \,=\,1$, but all other $M^p_n$ with $n\,\ge \, 2$=0, are equal to

$$\begin{aligned} M^p_n\left( Y, r\right) \,\,=\,\,\int ^{\epsilon \,+\,i\,\infty }_{\epsilon - i\,\infty }\frac{d \gamma }{2\,\pi \,i}e^{ \bar{\alpha }_S\chi \left( \gamma \right) \,Y}\,\Bigg ( e^{ \bar{\alpha }_S\chi \left( \gamma \right) \,Y}\, - \,1\Bigg )^{n - 1}\nonumber \\ \end{aligned}$$

(55)

which leads to the multiplicity distribution, which takes the form (see Eq. (41) and Ref. [73]):

$$\begin{aligned}&\tilde{P}_n\left( Y, r \right) \,\,\,=\,\,\int ^{\epsilon \,+\,i\,\infty }_{\epsilon - i\,\infty }\frac{d \gamma }{2\,\pi \,i}\, e^{-\, \bar{\alpha }_S\chi \left( \gamma \right) \,Y}\,\nonumber \\&\qquad \quad \qquad \quad \times \Bigg ( 1\, - \,e^{ -\, \bar{\alpha }_S\chi \left( \gamma \right) \,Y}\Bigg )^{n - 1} \end{aligned}$$

(56)

For $N\,\,=\,\,e^{ \bar{\alpha }_S\chi \left( \gamma \right) \,Y} \,\,\gg \,\,1$ we have

$$\begin{aligned}&\tilde{P}_n\left( Y, r \right) \nonumber \\&\quad =\,\,\int ^{\epsilon \,+\,i\,\infty }_{\epsilon - i\,\infty }\frac{d \gamma }{2\,\pi \,i}\,\exp \left( - z(\gamma )\,\,+\,\,\gamma \,\xi \right) \nonumber \\&\hbox {where}z\,\,=\,\,\frac{n}{e^{ \bar{\alpha }_S\chi \left( \gamma \right) \,Y} } \end{aligned}$$

(57)

Taking the integral over $\gamma$ using the method of steepest descent, using the diffusion approximation for the BFKL kernel (see Eq. (46)). The equation for $\gamma _\mathrm{SP}$ has the following form:

$$\begin{aligned}&2\,D\,Y z\left( \frac{1}{2}\right) \left( \gamma _\mathrm{SP} - \frac{1}{2}\right) +\,\xi \,=0~~~~\nonumber \\&\quad \hbox {with}~~~\left( \gamma _\mathrm{SP} - \frac{1}{2}\right) \,\,=\, - \frac{\xi }{2\,D\,Y\,z\left( \frac{1}{2}\right) } \end{aligned}$$

(58)

and the integral over $\gamma$ is

$$\begin{aligned}&\tilde{P}_n\left( Y, r \right) \,\,\,=\,\,\, \sqrt{\frac{\pi }{2\,D\,z\left( \frac{1}{2}\right) \,Y}}\, \frac{1}{ N\left( Y\right) }\,e^{ - z\left( \frac{1}{2}\right) }\, ~~~\hbox {with} ~~~~z\nonumber \\&\quad =\,\frac{n}{N\left( Y\right) }~~~~\hbox {and}~~~N\left( Y\right) \,\,=\,\,e^{\omega _0\,Y} \end{aligned}$$

(59)

considering $\xi ^2\Big /2\,D\,z\left( \frac{1}{2}\right) Y \,\ll \,1$. After normalization, we obtain that

$$\begin{aligned}&\frac{< n > \,\sigma _n}{\sigma _\mathrm{in}}\,\,\,=\Psi ^\mathrm{KNO} \left( z\right) \,\,\,=\,\,\, \sqrt{\frac{1}{\pi \,z}}\, \,e^{ - z\left( \frac{1}{2}\right) }\, \nonumber \\&~~~\hbox {with} ~~~~z\,=\,\frac{n}{N\left( Y\right) }~~~~\hbox {and}~~~N\left( Y\right) \,\,=\,\,e^{\omega _0\,Y} \end{aligned}$$

(60)

where $\Psi ^\mathrm{KNO}$ denotes the KNO function (see Ref. [89]).

It is worthwhile mentioning that the multiplicity distribution of Eq. (60) is different from Eq. (27) and

$$\begin{aligned} R\,=\,\frac{P_n^\mathrm{QCD}\left( Eq.~(60)\right) }{P_n\left( Eq.~(27)\right) }\,\,=\,\,\sqrt{\frac{1}{\pi \,z}} \end{aligned}$$

(61)

In Fig. 6 we compare the ALICE data [90] on multiplicity distribution with Eq. (59) and with Eq. (27). On can see that the agreement is good, and the difference between the above equations can be seen at large n. In describing the experimental data we use Eq. (59), which is derived at large z, for $z\,\ge \,3$. It worthwhile mentioning that the data of CMS [91] we have discussed in our paper [73].

Fig. 6

Multiplicity distribution of the charged hadrons in the central rapidity region. The solid line is the distribution of Eq. (59), while the dotted curve corresponds to Eq. (27). The data and the value of $N=12$ are taken from Ref. [90]

4.2 $P^{(2)}_n$ distribution for quarkonia production

The $P^{(1)}_n$ distribution , which we have discussed in the previous section, can be derived, using the double Laplace transform representation, both for $P^{(1)}_n\left( Y; r\right)$ and for $M^{(1)}_n\left( Y; r\right)$:

$$\begin{aligned} M^{(1)}_n\left( Y; r\right)= & {} \int ^{\epsilon \,+\,i \infty }_{\epsilon - i\,\infty } \frac{d \omega }{2\,\pi \,i}\nonumber \\&\int ^{\epsilon \,+\,i \infty }_{\epsilon - i\,\infty }\frac{d \gamma }{2\,\pi \,i}\,e^{\omega \,Y\,\,+\,\,\gamma \,\xi }\,\, m^{(1)}_n\left( \omega , \gamma \right) \nonumber \\ \end{aligned}$$

(62)

where $\xi \,=\,\ln \left( r^2 \Lambda ^2_\mathrm{QCD}\right)$.

Equation (54) in the $\omega$-representation, has the following form:

$$\begin{aligned}&\omega \,m^{(1)}_n\left( \omega , \gamma \right) \,\,\,=\,\,n\,\bar{\alpha }_S\chi \left( \gamma \right) \,m^{(1)}_n\left( \omega , \gamma \right) \,\,\nonumber \\&+\,\,\left( n - 1 \right) \,\bar{\alpha }_S\chi \left( \gamma \right) \,\,m^{(1)}_{n-1}\left( \omega , \gamma \right) \end{aligned}$$

(63)

with the solution:

$$\begin{aligned} m^{(1)}_n\left( \omega , \gamma \right) \,\,\,=\,\,\left( n\, -\,1\right) ! \prod ^{n}_{m=1} \frac{ \bar{\alpha }_S\chi \left( \gamma \right) }{ \omega \, - \,m\,\bar{\alpha }_S\chi \left( \gamma \right) } \end{aligned}$$

(64)

This solution gives $M^p_1\left( Y=0,r\right) \,\,\ne \,\,0$, while $M^p_n\left( Y{=}0,\right. \left. r\right) \,\,=\,\,0$ for $\,n\,\ge \,2$. The inverse Laplace transform leads to Eq. (55) for $M^p_n\left( Y, r\right)$ and Eq. (56) for $P^{(1)}_n\left( Y, r\right)$.

To find $P^{(2)}_n$ distribution we need to take into account that at $Y=0$:

$$\begin{aligned}&M^p_2\left( Y = 0, r\right) \,\,=\,\,1;~~~~~~~~M^p_1\left( Y = 0, r\right) \,\,=\,\,0;\nonumber \\&M^p_n\left( Y = 0, r\right) \,\,=\,\,0~~\hbox {for}~~n\,\,\ge \,\,3 \end{aligned}$$

(65)

One can see that the following: $m^{(1)}_1\left( \omega , \gamma \right)$ satisfies these conditions:

$$\begin{aligned} m^{(2)}_n\left( \omega , \gamma \right) \,\,\,=\,\,\left( n\, -\,1\right) ! \prod ^{n}_{m=2} \frac{ \bar{\alpha }_S\chi \left( \gamma \right) }{ \omega \, - \,m\,\bar{\alpha }_S\chi \left( \gamma \right) } \end{aligned}$$

(66)

The inverse Laplace transform with respect to $\omega$ leads to

$$\begin{aligned} M^{(1)}_n\left( Y; r\right)= & {} \int ^{\epsilon \,+\,i \infty }_{\epsilon - i\,\infty }\frac{d \gamma }{2\,\pi \,i}\,e^{\gamma \,\xi }\,\,\frac{1}{\gamma }\,\left( n - 1\right) \, e^{2\,\bar{\alpha }_S\,\chi \left( \gamma \right) \,Y}\nonumber \\&\times \Big (e^{2\,\bar{\alpha }_S\,\chi \left( \gamma \right) \,Y}\, - \,1\Big )^{n - 2} \end{aligned}$$

(67)

which gives the initial conditions of Eq. (65).

For $P^{(2)}_n$ we obtain

$$\begin{aligned} P^{(2)}_n\left( Y; r\right)&{=}&\int ^{\epsilon \,{+}\,i \infty }_{\epsilon {-} i\,\infty }\frac{d \gamma }{2\,\pi \,i}\nonumber \\&\times e^{\gamma \,\xi }\,\,\frac{1}{\gamma }\,\left( n - 1\right) \, e^{-\,2\,\bar{\alpha }_S\,\chi \left( \gamma \right) \,Y}\,\Big (1\, {-} \,e^{-2\,\bar{\alpha }_S\,\chi \left( \gamma \right) \,Y}\Big )^{n - 2}\nonumber \\ \end{aligned}$$

(68)

with $P^{(2)}_2\left( Y=0,r\right) \,\,=\,\,1$ and $P^{(2)}_n\left( Y=0,r\right) \,\,=\,\,0$ for $n\,\ne \,2$ at $Y=0$.

Repeating the same estimates as in Eq. (57), we obtain the KNO function,

$$\begin{aligned} \Psi ^{KNO}\left( z\right) \,\,=\,\,2\,\sqrt{ \frac{z}{\pi }}\,e^{-\,z} \end{aligned}$$

(69)

with the normalization $\int d z\,\Psi ^{KNO}\left( z\right) \,\,=\,\,1$.

Note that the ratio $\frac{ P^{(2)}_n\left( Y, r\right) }{P^{(1)}_n\left( Y , r\right) } \,\,=\,\,z$ for large z, as in Eq. (39).

5 Comparison with experimental data

In both ALICE [52,53,54,55,56] and STAR  [57, 58] experiments the following ratio is measured:

$$\begin{aligned} \frac{n_{J/\Psi }}{\langle n_{J/\Psi }\rangle } \,\,=\,\,F\left( \frac{n}{N}\right) \end{aligned}$$

(70)

where $N\,=\langle n \rangle$ is the average number of charged hadrons in the fixed rapidity window, and $\langle n^{J/\Psi }\rangle$ the average number of $J/\Psi$ which are measured generally speaking in a different rapidity window. It turns out that $F\left( \frac{n}{N}\right) \,\,\ne \,\,\frac{n}{N}$, but it is close to this when the rapidity windows are different. When both rapidity windows are the same, $F\left( \frac{n}{N}\right)$ shows much steeper dependence than $\frac{n}{N}$. In Ref. [60] the $J/\Psi$ production is considered as being proportional to the number of collisions, since it comes from short distances, while the production of hadrons is proportional to the number of participants (see Refs. [77,78,79,80]). However, in the framework of the CGC approach, the $J/\Psi$ production at high energies is proportional to the number of participants [61, 66,67,68,69] as it can be seen from Fig. 1.

The main ingredients for describing the experimental data are Eqs. (56)–(60) and Eqs. (68)–(69), as well as Eq. (37). Using these equation we can re-write Eq. (35) in the form:

$$\begin{aligned}&\frac{\frac{d \sigma _n^\mathrm{J/\Psi }}{ d y}}{ \frac{d \sigma ^\mathrm{prod. gl.}_n}{ d y}} \nonumber \\&\quad =\,\, \frac{\frac{d \sigma _\mathrm{incl}^\mathrm{J/\Psi }}{ d y}}{ \frac{d \sigma ^\mathrm{prod. gl.}_\mathrm{incl}}{ d y}} \frac{1}{1\,\,+\,\,\kappa }\Bigg ( \frac{n}{\langle n^{(1)}\rangle } \ \,\,\,+\,\,\, \kappa \frac{n}{\langle n^{(2)}\rangle } \, \, \frac{ P^{(2)}_n\left( N\right) }{ P^{(1)}_n\left( N\right) }\Bigg ) \nonumber \\&\quad =\,\,\frac{\frac{d \sigma _\mathrm{incl}^\mathrm{J/\Psi }}{ d y}}{ \frac{d \sigma ^\mathrm{prod. gl.}_\mathrm{incl}}{ d y}} \frac{1}{1\,\,+\,\,\kappa }\Bigg ( \frac{n}{N} \ \,\,\,+\,\,\, \kappa \frac{n}{2\,N } \, \, \frac{ P^{(2)}_n\left( N\right) }{ P^{(1)}_n\left( N\right) }\Bigg ) \end{aligned}$$

(71)

Using Eq. (71) we can find the experimental observable (see Refs. [53, 63])

$$\begin{aligned}&\displaystyle {\frac{dN_{J/\psi }/dy}{\langle dN_{J/\psi }/dy\rangle }\,\,=\frac{w\left( N_{J/\psi }\right) }{\left\langle w\left( N_{J/\psi }\right) \right\rangle }\,\frac{\left\langle w\left( N_{\mathrm{ch}}\right) \right\rangle }{w\left( N_{\mathrm{ch}}\right) }} \nonumber \\&\quad = \frac{d\sigma _{J/\psi }\left( y,\,\eta ,\,\sqrt{s},\,n\right) /dy}{d\sigma _{J/\psi }\left( y,\,\eta ,\,\sqrt{s},\,n=N\right) /dy}\nonumber \\&\quad / \frac{d\sigma _{\mathrm{ch}}\left( \eta ,\,\sqrt{s},\,n\right) /d\eta }{d\sigma _{\mathrm{ch}}\left( \eta ,\,\sqrt{s},\,n = N\right) /d\eta }\nonumber \\&\quad = \frac{\Bigg ( \frac{n}{N} \ \,\,\,+\,\,\, \kappa \frac{n}{2\,N} \, \, \frac{ P^{(2)}_n\left( N\right) }{ P^{(1)}_n\left( N\right) }\Bigg )}{ \Bigg (1\ \,\,\,+\,\,\, \frac{1}{2}\kappa \, \, \frac{ P^{(2)}_N\left( N\right) }{ P^{(1)}_N\left( N\right) }\Bigg )} \displaystyle { \xrightarrow {z = \frac{n}{N} \gg 1} \frac{ z\,\,+\,\,\frac{\kappa }{4}\,z^2}{1+\frac{\kappa }{4}}} \end{aligned}$$

(72)

In Fig. 7 we compare the experimental data with Eq. (72).

Fig. 7

Comparison Eq. (72) with the experimental data of the ALICE collaboration [52,53,54,55,56]. The solid line is the estimate of Eq. (72), and the dotted line is the linear dependence which stems from the contribution of a. a Shows the production of $J/\Psi$ in central rapidity region, while in b the estimates are shown for Eq. (72) with $\kappa$ calculated in leading order of perturbative QCD, and with $\bar{\alpha }_S= 0.15$

Fig. 8

The production of quarkonia from n parton cascades

One can see that this simple formula provides a fairly good description of the experimental data, for central production where $\kappa = 1$. However, the experimental data for forward production of $J/\Psi$ [52,53,54,55,56,57,58] show almost a linear dependance: $\frac{n_{J/\Psi }}{\langle n_{J/\Psi }\rangle }\,\,=\,\,z$. Indeed, in Eq. (72) the quadratic term is suppressed, since it is proportional to the value of $\kappa$, which is equal to (see Fig. 4 and Ref. [61] for the estimates).

$$\begin{aligned} \kappa \,\,\,\,=\,\,\,\left( \frac{Q^2_s\left( Y - y\right) }{Q^2_s\left( y\right) }\right) ^{ \bar{\gamma }}\,\,=\,\,e^{-\,2\, \bar{\gamma } \lambda \,y^*} \end{aligned}$$

(73)

where $y^*$ is the rapidity of the produced quarkonia in c.m.f. In leading order of perturbative QCD , in which we made all our previous estimates [1], $\bar{\gamma }\,\,=\,\,0.63$ and $\lambda \,=\,\bar{\alpha }_S\,\frac{\chi \left( \bar{\gamma }\right) }{\bar{\gamma }}\,\,\approx \,4.8\,\bar{\alpha }_S$. As is shown in Fig. 7b, the estimate in leading order describes the data quite well. However, we need to remember that the NLO corrections both to $\bar{\gamma }$ and to $\lambda$ are large.

6 Conclusions

In this paper we re-visited the problem of multiplicity distributions in high energy QCD, which we have discussed in Ref. [73] and found the distribution of Eq. (59). This distribution provides a better description of the experimental data at large multiplicities n, than Eq. (27), which has been discussed previously. We also suggest a different approach to the multiplicity dependence of quarkonia production. It should be stressed that our approach is based on the three gluons fusion mechanism of Fig. 1, and it differs from the description of Refs. [62, 63], since we did not assume the multiplicity dependence of the saturation scale. In our approach we assume, that the production of $J/\Psi$, which occurs at rapidity y, and the central production of charged hadrons, stem from the production of the same n-parton cascades, which are pictured in Fig. 8 as the production of n-gluon ladders. Solving the QCD cascade equation, we found the multiplicity distribution both for the cascade of Fig. 8a (see Eq. (68)–(69)) and for the cascade of Fig. 8b (see Eq. (56)–(60)) .

In Fig. 8 one can see that $J/\Psi$ can be produced from each of n-ladders, leading to the cross section, which is proportional to n (see Fig. 8a). This mechanism is shown in Fig. 3a. However, $J/\Psi$ can be created from merging of two ladders (see Fig. 8b) , which gives a cross section $\propto \,\,n^2$, and corresponds to Fig. 3b. Note, that the production of the hadrons in both cases can be found from Eq. (37). Taking into account that the average number of gluons (hadrons) for the mechanism of Fig. 8b is two time larger than for Fig. 8a, we infer that Fig. 8 leads to the simple Eq. (72).

It should be stressed that this equation is heavily dependent on the three gluon fusion mechanism, but does not depend on the details of the cross section of quakonia production. In particular, as we have mentioned above, we do not use the dependence of the saturation scale on the multiplicity of the produced gluons. This means that the non-linear dependence of $J/\Psi$-production on multiplicity of charged hadrons, can stem from sources other than the dependence of the cross section on the saturation scale. Actually, this statement follows directly from the fact that 1 + 1 RFT generates the non-linear dependence on n.

It should be noted that if we assume, that in addition to the three gluon fusion mechanism we have the production of $J/\Psi$ from color-singlet model, we will not obtain agreement with the description of the experimental data in Fig. 7. In particular, in the colour-singlet model(CSM), which is able to describe the experimental data on $J/\Psi$ inclusive production [92] we expect that $\frac{dN_{J/\psi }/dy}{\langle dN_{J/\psi }/dy\rangle }\,\propto \,n$. In our approach, any non-linear dependence in the CSM is an indication that we have to develop a different approach to QCD cascade or suggest another model for production of soft hadrons.

As an aside, we note in [62, 63, 93] it was assumed that the $J/\Psi$-production results from the inclusive diagrams of Fig. 3a and b. This is erroneous, as at arbitrary n it is necessary to include the production of many partonic showers as illustrated in Fig. 4a and b. The difference to the percolation approach [60], lies in our hypothesis that both the production of $J/\Psi$ and the charged pion stem from short distances of the order of $r \, \propto \,1/Q_s$, and are determined by physics controlled by the CGC effective theory. The gluon jets with transverse momentum $Q_s$, decay into charged pions (see Ref. [63] for details). The non-linear dependence of production of $J/\Psi$ is due to the three Pomeron fusion mechanism.

Following the advice of our referee, we add two comments. First, the open charm production is enhanced at large multiplicities [94]. In spite of the fact, that three gluon fusion gives a sizable contribution to inclusive production  [93], the two gluon fusion leads to main contribution to this process. We are planning to discuss the multiplicity distribution in this process in our future publication. The discussion in Ref. [93] cannot be trusted since it does take into account the production of all n ladders in Fig. 8. The failure of our approach to explain the experimental data will indicate, that we need to introduce the dependence of the saturation scale on the multiplicity of produces hadrons (see Refs. [63, 82, 93]). Second, the production of $\Upsilon$-meson stems from the short distances, much shorter than $1/Q_s$. The three gluon fusion gives the main contribution to $\Upsilon$ inclusive at high energies, but the corrections due to CSM at the LHC could be sizable. We plan to consider this process in our future papers.

In spite of the good description of the experimental data for the quarkonia production integrated over the transverse momenta ($p_T$), we cannot explain at present, why the data at fixed $p_T$  [52], shows a steeper dependence on n than the integrated data. Certainly, this problem will be the main subject of our further attempts to understand the multiplicity dependence of quarkonia production.

Concluding, we wish to stress that we are only in the beginning of the theoretical understanding of multiplicity distributions, and this paper is the first attempt to give the self consistent description of this distribution in the framework of perturbative QCD with the simplified assumption about confinement quarks and gluons: the local hadron–parton duality.