Four-pomeron vertex

The four-pomeron vertex is studied in the perturbative QCD. Its dominating terms of the leading (zeroth and first) orders in the coupling constant and subdominant in the number of colors are constructed. The vertex consists of two terms, one with a derivative in rapidity $\pd_y$ and the other with the BFKL interaction between pomerons. The corresponding part of the action and equations of motion are found. The iterative solution of the latter is possible only for rapidities smaller than 2 and quite large coupling constant $\alpha_s$, of the order or greater than unity, when the quadruple pomeron interaction is relatively small. Also iteration of the part with $\pd_y$ is unstable in the infrared region and compels to introduce an infrared cut. The variational approach with simple trying functions allows to find the minimum of the action at $\alpha_s$ of the order 0.2 and rapidities up to 25. Numerical estimates for O-O collisions show that actually the influence of the quadruple pomeron interaction turns out to be rather small.


Introduction
For many years the high-energy behavior in the QCD has been the subject of intensive study both in the context of the so-called JIMWLK approach (see e.g. [1] and references therein) and of the BFKL-Bartels approach based on summation of the diagrams constructed for the interaction of reggeized gluons ('reggeons') and summarized in the effective QCD reggeon action ( [2]). For the interaction of a point projectile with a large nucleus and in the approximation of a large number of colors N c in both approaches one arrives at a simple closed equation, the Balitski-Kovchegov (BK) equation [3,4], which actually sums the fan pomeron diagrams constructed with the BFKL Green functions and the triple pomeron vertex. In our papers [5,6,7] this equation has been generalized for collisions of two heavy nuclei. Unlike the BK equation the analogous equation for AB collisions is no more an equation for evolution in rapidity but rather the one with boundary conditions at initial and final rapidities and so much more difficult, as illustrated by quite few attempts at its solution [8,9,10].
Both the BK equation and its generalization made so far are based on the pomeron interaction via the triple pomeron vertex. While for the BK equation in the adopted approximation (lowest order, absence of loops) it is sufficient, it is not so for the interaction of heavy nuclei where also interactions mediated by the quadruple-pomeron vertex may be important. Its appearance can be traced to the form of the gluon production in [11], which obviously included production from the quadruple pomeron vertex. Studies in the drastically simplified ("toy") models without transverse dimensions have shown a strong influence of the quadruple pomeron interactions on the high-energy behavior [12,13]. For this reason it is interesting and important to study the quadruple pomeron vertex in the QCD, which is the subject of this article. In the equations for the nucleus-nucleus amplitude it will appear as a new interaction.
We recall that these equations are obtained from the non-local action where all parts are obtained from the bilocal fields ψ(y; r 1 , r 2 ) and ψ † (y; r 1 , r 2 ) describing incoming and outgoing pomerons and depending on rapidity y and two spatial points r 1 and r 2 of the two reggeons of which they consist. In Eq. (1) part S 0 describes free pomerons Here z combines the two points r 1 and r 2 : z = {r 1 , r 2 }.
Hamiltonian H can be taken in the form symmetric in the incoming and outgoing pomerons. Then the pomeron P (k 1 , k 2 ) ("semi-amputated") in the momentum space is related to the standard BFKL pomeron P BF KL (k 1 , k 2 ) as P (k 1 , k 2 ) = k 1 k 2 P BF KL (k 1 , k 2 ) and where ω(k) is the reggeized gluon trajectory minus one. The free action defines the pomeron propagator G(y; z, z ′ ) which satisfies (∂ y + H)G(y; z, z ′ ) = δ(y)δ 4 (z − z ′ ).
Part S E determines the interaction of the pomerons with the participant nuclei. It is assumed that the sources for the projectile B and target A are different from zero at rapidities y = Y B and y = Y A respectively: If Y B > Y A then ψ(y, z) = 0 at y < Y A and ψ † (y, z) = 0 at y > Y B . So in the action the integration over y actually extends over the interval Y A < y < Y B . The part S 3 describes the interaction between the pomerons. Both in the BK equation and its generalization for AB scattering made so far it was taken as the triple pomeron interaction. Then in the configuration space it has the form This expression looks quite complicated due to nonlocal operators T . However in diagrams it is coupled to pomeron propagators, which also contain these operators so that in the end all the non-locality is eliminated.
Our aim is to find an additional part for interaction S 4 which comes from the quadruple pomeron interaction. This part contains various contributions with different orders in the coupling g and number of colors N c . So our first task will be to study these orders, which will be done in the next section. The third section will be devoted to the derivation of the leading contribution for S 4 . In the fourth section we shall discuss the changes in the coupled pomeron equations due to quadruple interaction in the no-loop approximation. In the fifth section we shall estimate the actual contribution of this interaction to the total O-O cross-section. Some conclusion will be presented in the last section.

Orders of magnitude
To estimate the orders of magnitude in the pomeron interactions one can use the simplest diagrams with projectile and targets represented by simple quarks attached to the lowest order pomerons (just two-gluon exchange). Additional reggeon interactions in the BFKL approach will contribute contributions of the relative orders (α s N c Y ) n ∼ 1 where Y is the overall rapidity. One should take into account that diagrams with a vertex, that is with an interaction between reggeons at a fixed rapidity y, will contain factor Y due to integration over y. In the pure perturbative approach one should additionally take into account the coupling to external quarks. So the rule is to take lowest order diagrams with quarks as participants and divide by factor (α s N c ) n P where n P is the number of participants. In this way for the triple pomeron vertex we find the total order α 5 s N 4 c and for the vertex, dividing by (α s N c ) 3 , we obtain the correct factor α 2 s N c . The quadruple pomeron vertex, unlike the triple one, does not change the number of reggeons. The diagrams of the lowest order include no pomeron interaction at all. The dominant diagram is of course the disconnected one Fig. 1,A of the total order α 4 s N 4 c . With N P = 4 this gives its final order unity.
There is another diagram with no pomeron interaction but with a redistribution of color in the collision Fig. 1,B. It will give a contribution to the quadruple pomeron interaction. Its order is α 4 s N 2 c . This gives order 1/N 2 c to the resulting quadruple pomeron vertex. Next we consider diagrams with reggeon interactions. As mentioned they will contain factor Y . The lowest order diagram is shown in Fig. 1,C. Its order is α 5 s N 3 c Y . Dividing by (α s N c ) 4 we obtain the factor for the corresponding four-pomeron vertex α s /N c . Factor Y ∼ 1/α s N c converts the total contribution into 1/N 2 c , so that effectively this diagram gives a contribution of the same order as Fig. 1 Diagrams with transition of 2 reggeons into 3 Fig. 1,D and 3 reggeons into 3 Fig. 1 which follows from α s N c y ∼ 1 at Y >> 1 assumed in the BFKL approach. So to conclude, the dominant contribution to the quadruple interaction comes from diagrams of the type Fig. 1,B and C. Vertices corresponding to diagrams Fig. 1 The BFKL approach assumes α s N C << 1 but we take N c >> 1, so this relative order is in fact undetermined.
3 Quadruple pomeron vertex in the lowest approximation 3.1 Fig. 1

.B
As follows from the previous section the dominant contributions to the 4-pomeron interaction come from the diagrams of the type shown in Figs, 1,B and C. We start from the first diagram Fig. 1.B. Taking for the participants two incoming and two outgoing pomerons we shall actually consider the diagram shown in Fig. 2. To write down the expression for this diagram it is convenient to make use of the multirapidity formalism for the reggeon dynamics introduced in [14]. In this formalism resembling the old Gribov technique, each reggeon is assumed to possess its own "energy" ǫ apart from its transverse momentum k with the propagator This energy is assumed to be conjugated to the rapidity variable y, which is established by the transition of time t to rapidity y according to it → y. At each interaction both energies and transverse momenta are conserved. In this formalism the pomeron with energies and momenta of its two component reggeons ǫ 1 , k 1 and ǫ 2 , k 2 can be presented as where E 12 + ǫ 1 + ǫ 2 is the total pomeron energy. The "normal" pomeron P (E 12 , k 1 , k 2 ) depending on its own energy E 12 conjugated to its rapidity is given by .
Now we turn to the diagram in Fig.2. We assume that the two incoming pomerons have their energies E 12 and E 34 and the two outgoing pomerons E 13 and E 24 . The conservation law gives factor 2πδ( which will be suppressed in the following. Conservation of energy at each pomeron further restricts the reggeon energies to a single energy, say, ǫ 1 with the rest determined as follows The diagram also depends on fixed total momenta of the pomerons k 12 = k 1 + k 2 and so on with k 12 + k 34 = k 13 + k 24 . As a result there is only one independent transverse momentum, say k 1 and the rest given as So the diagram will be given by the integral over energy ǫ 1 and transverse momentum k 1 . Suppressing the latter integration we find Here F 12 = F (E 12 ; 1, 2), for brevity we denote the transverse momenta by just their numbers: 1 = k 1 and so on. It is assumed that each F carries factor α c N c which produces the overall factor 1/N 2 c . All ω's are assumed to have a small positive imaginary part. Calculation of the integral over ǫ 1 gives In terms of pomerons P (E; k 1 , k 2 ) the diagram will be given after the transverse momentum integration of where P 12 ≡ P (E 12 ; 1, 2) and so on. Passing to dependence on rapidities we have D 1 (y 12 , y 34 , y 13 , y 24 ) = 1 N 2 where we discriminate between the incoming and outgoing pomerons, P 12 , P 34 and P 13 , P 24 respectively. We introduce the vertex rapidity y = it by presenting and for each pomeron use with the subsequent passage t → −iy, in which an extra factor (−i) appears, see [14]). In this way we ultimately get D B (y 12 , y 34 , y 13 , y 24 ) = 1 N 2 × dyP (y 12 − y; 1, 2)P (y 34 − y; 3, 4)P (y − y 13 ; 1, 3)P (y − y 24 ; 2, 4).
Note that P (y) has a jump at y = 0. One can put the derivatives inside the integral and rewrite (21) in the form So the contribution from the diagram Fig. 1,B splits into two parts ω i dyP (y 12 − y; 1, 2)P (y 34 − y; 3, 4)P (y − y 13 ; 1, 3)P 24 (y − y 24 ; 2, 4). (24) The term D B is infrared divergent due to trajectories ω i . It has order α s N c and is accompanied by factor Y in the diagram, so it has the same order as the diagram in Fig. 1,C to be considered later. In the sum with the latter the infrared divergence will be canceled.
The term D (1) 1 is infrared finite. To see its order it is convenient to transform it using the BFKL equation (5). The incoming pomeron can be presented via the Green function as where K combines k 1 and k 2 , ρ 12 is the color source and in the second equality ρ and G(y) are considered as operators in the K space, conjugate to the coordinate z-space introduced in the Introduction. The Green function satisfies Eq. (5), which is the operatorial notation is From this equation we find where the momentum or coordinates are suppressed and it is assumed that Hamiltonian H 12 acts on them. Similarly for the second incoming pomeron For the outgoing pomeron we have so that we get in the same way and ∂ y 24 P 24 (y − y 24 ) = −ρ 34 δ(y − y 24 ) + H 24 P 24 (y − y 24 ).
As a result the contribution D 1 in its turn splits into two parts. The first part comes from the differentiation of δ functions in rapidities: Graphically it corresponds to Fig. 3.
Restoring the transverse momentum integration and introducing independent momenta for the outgoing pomerons we finally have and The diagram shown in Fig. 1,C includes interaction between reggeons 2 and 3, which with the external pomerons is shown in Fig. 4. Its contribution should be summed with a similar contribution with the gluon interaction between reggeons 1 and 4. As before we denote initial total energies as E 12 and E 34 and the final energies as E 13 and E 24 . We again suppress the energy conservation factor 2πδ(E 12 + E 34 − E 13 − E 24 ). The diagram contains two loops and so energy integrations can be taken over ǫ 1 and ǫ 4 . Energies of the gluons 1, 2, 3 and 4 before the interaction are After the interaction gluons 2 and 3 have their energies E 24 −ǫ 4 and E 13 −ǫ 1 respectively and change their momenta k 2 → k ′ 2 and k 3 → k ′ 3 . Note that the interaction connects different color configurations with gluons from the projectile forming colorless pairs (1,2) and (3,4) and from the target forming colorless pairs (1,3) and (2,4). As a result the interaction has the opposite sign as compared to the one inside the pomeron [14]. If the pomeron Hamiltonian is then the interaction in Fig. 4 is just −V . So in terms of amputated pomerons F the contribution from Fig. 4 before the transverse integrations is .
Integration over energies factorizes into two integrals: and a similar integral over ǫ 4 which gives .
So passing to the pomerons we get The second amplitudeD C with the interaction between the gluons 1 and 4 is Both D C andD C are infrared divergent. Using the fact that the BFKL Hamiltonian is infrared stable as a whole, we can separate the divergence into the reggeon trajectories, presenting (see (37)) Note that for arbitrary momenta of the outgoing pomerons in the transversal integral D C has to be multiplied by (2π In both cases terms with ω i acquire the product of all 4 delta-functions between initial and final momenta. Then one finds that the terms with delta-functions cancel with the similar terms in D After this cancelation the resulting contribution from the diagrams with one reggeon interaction can be written as Transition to rapidities, as before, will give D 1 (y 12 , y 34 , y 13 , y 14 ) = − 1 N 2 c dyP (y 12 − y; 1, 2)P (y 34 − y; 3, 4)P (y − y 13 ; 1 ′ , 3 ′ )P (y − y 24 ; 2 ′ , 4 ′ )

The quadruple pomeron vertex and action
By definition the quadruple pomeron vertex Γ corresponds to a contribution generated by the coupling of four pomeron propagators, that is the Green function G, at rapidity y of the structure illustrated in Fig. 6. In the vertex pairs of momenta are grouped to correspond to colorless states. We also define it with the extracted δ function corresponding to conservation f the total momentum. Inspecting our results in the previous section we can identify contributions D B and D 1 with this structure. In the sum they give the vertex Figure 6: Four pomerons coupled with the quadruple vertex In accordance with (47) the effective action of pomeron fields acquires new interaction terms (note that the derivative acts on all field variables in (48)). and One can write down the quasi-classical equations following from δS/δψ = δS/δψ † = 0. To do this in the coordinate space one has to Fourier transform ψ and ψ † in (60) taking into account also the corresponding transformation of H 23 and H 14 . However these general equations have little sense since in applications to the nucleus-nucleus collisions the pomerons have zero total transverse momenta. In this case the equations substantially simplify (see [5,8]) and will be studied in the next section. The general, nonforward form of (60), whether in the coordinate or momentum space, is in fact essential only for the study of conformal invariance and calculating loops, the latter task remaining still very remote in future.

Forward case
In applications to the nucleus-nucleus scattering and in absence of loops all interacting pomerons have their total transverse momentum equal to zero. So the pomeron wave functions depend on the single transverse variable r or k in the coordinate or momentum spaces.
In terms of these pomeron fields the action with only the triple-pomeron interaction was introduced in [5,8].
The free action is then where H is the forward Hamiltonian (4).
Momenta k and k ′ are two-dimensional transversal and K is a differential operator in k Appearance of the operator K in (51) as compared to (2) is due to the transition of the semi-amputated pomerons in the coordinate space to the momentum space. The interaction part of the action describes splitting and merging of pomerons: The external action is where w A,B as before describe the interaction of the pomerons with the projectile and target. At fixed transverse point b A in the nucleus A coupling w A (k) is supposed to be proportional to the profile nuclear function T A (b A ) and its dependence on k is determined by the gluon distribution in the nucleon. Fields φ and φ † as well as w are dimensionless. At fixed b A and b B action S is dimensionful: dim S(b A , b B ) = 2. It becomes dimensionless after integration over b A at fixed impact parameter b, so that Note that in (51) it is assumed that the integration in y goes over the whole axis: −∞ < y < +∞ with the conditions that φ(y, k) = 0 at y < y min + 0 and φ † (y, k) = 0 at y > y m ax = Y . At y = 0 and y = Y the piece with ∂ y contributes terms proportional to δ(y) and δ(y − Y ) due to jumps of φ and φ † .
Turning to the quadruple pomeron interaction we have the same diagrams shown in Fig.  1 B and C which after the cancelation of the terms with ω give the contributions (34) and (44). Taking y 12 = y 34 = Y and y 13 = y 24 = 0 we find for the forward direction and For clarity we denote the upper and lower pomerons in our figures by subindexes u and l Note that in (57) H is not the forward Hamiltonian: From these expressions we can read the corresponding quadruple pomeron vertex and action S 4 . From D B we find the vertex Γ with the contribution to the action (60) For the equations of motion, which for φ(y, k) is obtained after differentiation δ/δφ † , one gets a contribution δ δφ † S The second contribution to the action from the quadruple pomeron interaction comes from (57) The classical equations of motion which follow, multiplied by (1/2)K −1 from the left, are and − ∂ ∂y From the δ-like dependence on y of the external sources it follows that the equations can be taken homogeneous in the interval 0 < y < Y , action of the external sources substituted by the boundary conditions in rapidities We noted in our previous papers that already with S 3 the equations of motion are no more pure evolution equations, since the initial conditions for φ(y) and φ † (y) had to be given at different rapidities, y = 0 for φ and y = Y for φ † . Inclusion of quadruple action S (1) 4 further complicates the situation, since according to (65) the initial conditions now each contain both φ and φ † .
To see the orders of magnitude of different terms in these equations it is instructive to rescale the pomeron fields φ and φ † , Hamiltonian and rapidity as whereᾱ = α s N c /π Then the equations take the form where and We observe that the on the right-hand side the free and triple interactions have the same order unity (actuallyᾱ in the original rapidity, which corresponds to the BK equation). The quadruple interaction has the orders 1/ᾱ 2 , which may be smaller or larger as compared with the two first terms depending on the relation between the small α s and large N c .

Cross-sections
At fixed overall impact parameter b and rapidity Y the nucleus-A-nucleus-B total crosssection is given by In the perturbative QCD the eikonal function T is given by the sum of all connected diagrams constructed of BFKL pomerons, which interact between themselves and with the participant nuclei, so that T is just the action S for fixed Y and b In the quasi-classical approximation (without loops) action S is to be calculated from the sum of (51), (54), 60), (62) and(55) with the solutions of evolution Eqs. (63) and (64), φ cl (y, k) and φ † cl (y, k). Using the equations of motion one can somewhat simplify the expression for S. Indeed multiplying the first equation by 2Kφ † , the second one by 2Kφ, integrating both over y and k and summing the results one obtains a relation which is valid for the classical action, that is, calculated with the solutions of Eqs. (63) and (64). Note that this condition is fulfilled for any form of interactions in S 2 , S 3 S 4 and S E , since it depends only on the number of field operators in them and hermiticity. Using this relation we can exclude, say, S 4 from (12) to find

Numerical estimates
As mentioned in our old paper [8] one can try to find the solution of our problem by two methods. One of them is to solve Eqs. (67) and (68) by iterations. In the simple case when by symmetry φ † (y, k) = φ(Y − y, k) in the first iteration one starts from Eq. (67) and solves it for φ(y) choosing for φ(Y − y) some simple initial value, typically just zero. Since in this case both the crossed term 2K −1 φ † (y, k)Kφ(y, k) and the term coming from S 4 are both equal to zero, the equation turns out to be the standard BK equation and its solution is just the sum of fan diagrams attached to the projectile. However at the next step one finds non-zero φ † (y) as φ(Y − y) from the obtained solution, so that both terms coming from S 3 and S 4 become present. Now on the right-hand side there appears a term which itself contains ∂ y φ, so that the equation is To solve it we divide the interval 0 < y < Y in steps ∆ and find at each step that is, using ∂ y φ obtained at the previous step on the right-hand side. One expects that with a sufficiently small ∆ the derivative ∂ y φ will not significantly change between steps and one will have good convergence. Solving in this way the equation one once again determines φ(y) and φ † (y) = φ(Y − y) and starts the third iteration of Eq. (67) and so on. If the iteration procedure is convergent one determines the solution φ(y, k) and φ † (y.k) = φ(Y − y, k) together with their derivatives in rapidity. Then one can find the external couplings w A and w B from (65). So in this approach these couplings are not fixed from the start but found after the equations for φ and φ † are solved. This circumstance has to be taken into account comparing these results with the pure perturbational calculations in Sec. 3. Unfortunately, as was mentioned in [8] this method has a very narrow region of convergence in A and Y , which shrinks with the growth of atomic number and rapidity. Without the quadruple interaction for O-O collisions it converges up to rescaled rapidityȳ = 2, which forᾱ = 0.2 corresponds to rather small actual rapidities less than 5. In attempt to consider higher rapidities and atomic number in [8] we recurred to another method trying to find the minimum of the action by a direct variational approach choosing some trying functions for φ and φ † . Even with very primitive trying functions with a varying ∆ in [8] we could find a minimum for the action for very large atomic numbers and rapidities, which lead to reasonable results for the eikonal function and cross-section. Comparison with the exact minimum found by the solution of the equation of motion at rapidities where such solution could be found by interactions showed that the precision of such variational approach was about 20%. In the present calculations our primary task was to see the influence of inclusion of the quadruple pomeron interaction. S 4 . We considered both methods, the solution of the quasiclassical equations Eqs. (67) and (68) by iterations and the direct variational approach. In both approaches we used the initial wave function which gave the best results in our old iterational solution of (67) and (68) where σ 0 = 20.8 mb and k is in GeV/c. This function is very similar in its behavior to the often used Golec-Biernat distribution but is simpler and leads to somewhat larger region of convergence. One can find some details of the calculations of different pieces of the action in Appendix.
In the iterative procedure to solve (67) and (68) we were confronted with a new problem related to the interaction containing the derivative ∂ y from S (1) 4 . Its contribution was found to lead to very high values of φ(y, k) in the infrared region, which prevented iterations already at very small values of y. Note that this problem does not exist for the second quadruple interaction coming from S (2) 4 which admits iterations at all k. To overcome this difficulty we had to introduce an infrared cut at the minimal value k min = exp t min for which iterations were possible. It turned out that t min ≃ −5. Even with this limitation the iteration procedure was found to be convergent practically in the same region as without S 4 only when the quadruple interaction was taken very small. Its relative magnitude is determined by factor 1/ᾱ 2 . For O-O collisions at the overall rapidity Y = 2 the iterative procedure turned out to be convergent only for quite large values of the coupling constant α ≥ 1. In our calculations we took the lowest possible valueᾱ = 1.1.
The gluon density at intermediate values of rapidity 0 < y < Y coming from both colliding nuclei can be determined as [8] dxG(x, k) A more illustrative comparison is done in Fig. 8. One observes first of all that at small momenta the contribution from the derivative coupling S blows up, which only exhibits the region of bad convergence. One cannot trust the obtained values of the gluon density in this boundary region. One rather expects the true density to behave roughly as with only S (2) 4 , the contribution from S (1) 4 only slightly lowering the density, as one can conclude from the figures at momenta immediately before the peak. At the physically important momenta in the region of the peak itself and higher the contribution from the quadruple interaction turned out to be very small due to the assumed high value ofᾱ.
The last conclusion is confirmed by the calculation of the eikonal function. In Fig 9 we present the eikonal at b = 0 for Y = 2, which is obtained from the full action (74) with the solution of equations (67) and (68)after integration over all b A = b B . The difference between the present and absent S 4 is smaller than 1% (and cannot be shown in the plot).
To conclude with the iterative solution of the quasi-classical equations for the fields, we kave found that it turns out to be dangerous in the infrared region. Probably the equations, being highly non-linear possess solutions of a different type, which avoid these difficulties but cannot be found by iterations.
Passing to direct variational search for the minimum of the action, we recall that without the quadruple interaction one always found the minimum using the trying function (77) with a variable ∆. Now we repeat this exercise with the quadruple interaction included. We drop the infrared cut and take a more reasonable value ofᾱ = 0.2, which of course substantially enhances the quadruple interaction. We find that inclusion of S 4 allows to find the minimum of the action up to rescaled rapidities 5 that is up to natural rapidities 25 but not to higher rapidities. So S 4 leads to certain restriction on the rapidity range, which however is physically not so bad. To simplify calculations we choose to calculate the action in the approximation of a constant T A (b A ) for each nucleus and afterwards obtained the eikonal function T(b) integrating the action in the overlap area of the collision depending on the impact parameter b. The obtained eikonal for central O-O collisions T (0) is illustrated in Fig 10 (lower curve). To compare we also present T (0) for the case without S 4 (upper

Conclusion
We have derived the quadruple pomeron interaction in the lowest order, which is essential in the interaction of nuclei between themselves. It consists of two terms, one with the derivative in rapidity and the other with the BFKL interaction between the pair of pomerons. Both are infrared safe. Their relative order to the triple pomeron interaction is 1/(α s N 2 ). So the quadruple interaction is much smaller for fixed α s and large number of colors. However in the opposite case with fixed N c and small coupling constant it may become comparable or even greater than the currently used triple interaction. Illustrative and very approximate calculations of the O-O cross-sections with the quadruple interaction together with the triple one and quite largeᾱ = 1 showed instability of the iterative procedure in the infrared region but a workable variational approach. We find a noticeable damping of the gluon density at small momenta, although at momenta greater than 1 GeV/c the change is quite small. In the variational approach one finds quite a substantial fall of the eikonal function at all impact parameters due to the quadruple interaction. More precise calculations of the AA collisions are needed to make a final conclusion of the importance of the quadruple interaction for the actual processes.

Appendix. Some details of the calculation
We use the logarithmic variable t putting k = exp(t) where k is measured in units provided by the initial condition φ(y, k) y=0 = φ 0 (k). Integration over k transforms as The integral was done by the Simpson numerical integration in n = 400 points. Without infrared cut satisfactory results were achieved for t min = −20, t max = +20. The infrared cut needed for iterations with S 4 was taken t min = −5.

Equation
The equation to solve can be rewritten in the form So after angular integration we get where z = t − t 1 , φ 1 = (∂/∂t)φ, φ 2 = (∂/∂t) 2 φ. The two quadruple interaction terms are where and Doing the angular integration and passing to integration variable t 1 we obtain, similarly to (82): where A is the same as in (83) and Finally the external term is T E (y, k) = 1 2 K −1 w A (k)δ(y).
To obtain our final expressions for calculation we first rescale φ, φ † and y as indicated in (66) and afterwards divide by 2α sᾱ . Then the coefficients in T 0 and T (1) 3 turn to unity and in T (2) 3 to 1/2. The coefficient in T 4 is divided by 4α 2 sᾱ . so we get for the total coefficient 1/(2π 2ᾱ2 ). Finally the rescaling of y produces factorᾱ in T E so that the coefficient becomes equal to α s .

Action
The action is given by Eqs. (51), (54), (60), (62) and (55) and is a sum Here we present it in the explicit form.
We begin with S 0 which is given by (51). It contains two terms: one with the derivative S  with H. We split operator K into the product K = L † L. In the logarithmic variable t one finds L = ∂ 2 /∂t 2 . So after the angular integration we get the part S where we again denoted φ 2 = ∂ 2 φ/∂t 2 . The integration over y is originally performed in the region −∞ < y < +∞ and so takes into account the jumps at y −0 and y = Y . Contribution from these jumps gives the boundary term where we have taken into account that φ † (0) = φ(Y ), The second term is S (2) 0 =ᾱ π dtdt 1 k 2 φ † 2 (y, k) B(k, k 1 )φ 2 (y, k 1 ) + A(k, k 1 )φ 2 (y, k) .
Next S (1) 4 = 1 πN 2 c dy dtk 2 φ † 2 φ∂ y φ. (105) Presence of the derivative ∂ y leads to the contribution from the jumps at y = 0 and y = Y and generates a new boundary term The second part is where w A is given by (99).