QCD unitarity constraints on Reggeon Field Theory

We point out that the s-channel unitarity of QCD imposes meaningful constraints on a possible form of the QCD Reggeon Field Theory. We show that neither the BFKL nor JIMWLK nor Braun’s Hamiltonian satisfy the said constraints. In a toy, zero transverse dimensional case we construct a model that satisfies the analogous constraint and show that at infinite energy it indeed tends to a “black disk limit” as opposed to the model with triple Pomeron vertex only, routinely used as a toy model in the literature.


Introduction
Reggeon Field Theory (RFT) is an effective theory for description of hadronic scattering in QCD at asymptotically high energies. The basic ideas of RFT go back to Gribov [1], and have been developed over the years in the context of QCD . In its modern form, the QCD RFT in a certain limit has been identified [35] with the so called JIMWLK evolution equation [36][37][38][39][40][41], or Color Glass Condensate (CGC) [42][43][44]. The relevant limit is when a perturbative dilute projectile scatters on a dense target.
Subsequently further relation between the CGC based approach and the RFT was explored. In particular recently we have shown that one can generalize the JIMWLK Hamiltonian consistently in the regime where large Pomeron Loops are important [45]. This regime includes the evolution of an initial dilute-dilute scattering to large rapidities,

JHEP08(2016)031
where at any given intermediate rapidity at most one of the evolved systems is dense. In this regime only Pomeron (Reggeon) splittings are important close to either one of the colliding objects, and one can write down a Hamiltonian, which encompasses both JIMWLK, and its dual (KLWMIJ [46]) evolution. The Hamiltonian in this regime contains the two triple Pomeron vertices, and (in the large N C limit) is the CGC equivalent of the Pomeron Lagrangian proposed by Braun [30][31][32] for description of the problem of scattering of large (but dilute) nuclei. So far, neither CGC nor RFT has been formulated in the most general case of scattering of two dense objects, although some work in this direction has been done [47][48][49][50][51][52].
There are some significant differences between the original Gribov RFT framework and its QCD incarnation. The original reggeons in Gribov's RFT are colorless, whereas the effective high energy degrees of freedom in QCD are frequently colored, such as reggeized gluons [53][54][55] or Wilson lines. It must be possible "to integrate over the color" and reformulate QCD RFT in terms of color neutral exchange amplitudes, such as BFKL Pomeron [2,3], however this has not yet been done explicitly. QCD RFT in addition to the Pomeron contains higher order colorless Reggeons, such as quadrupoles and higher multipoles. Whether these higher Reggeons significantly affect high energy behavior of QCD amplitudes is not known at present. Finally, even if the higher Reggeons can be discarded, it is not known whether the effective Pomeron Field Theory has a finite number of transition vertices. The large N C limit of high energy QCD is a convenient setup for the study of these questions. In this paper we stick to the large N c limit and in fact restrict ourselves even further by considering the dipole model approach [7][8][9], in which the Pomeron is the only relevant degree of freedom at high energy.
It is clear by now that the CGC formalism conceptually provides a direct route to derive the Reggeon Field Theory from the underlying QCD. Due to this direct connection, one expects that it should be possible to understand some general features of RFT that are required by QCD. The current paper is devoted to discussion of how the unitarity of QCD as a fundamental field theory exhibits itself in the RFT framework. To be precise we will be discussing the s-channel unitarity.
Our motivation to consider this question largely comes from earlier studies of Pomeron Lagrangian proposed by Braun [30][31][32] for scattering of large (but dilute) nuclei. The Lagrangian incorporates the BFKL dynamics in the linear regime and contains two symmetric triple Pomeron vertices. When the scattering amplitude is evolved within this framework to high enough rapidity, it exhibits paradoxical behavior: classical solutions to the equations of motion bifurcate beyond some critical rapidity Y c [56] and the dependence of the Pomeron amplitude on rapidity becomes unphysical. One is then left to wonder whether this peculiarity is a consequence of a possible non-unitarity of the Braun evolution.
The aim of this paper is to formulate the requirements of QCD (s-channel) unitarity in the (Pomeron) RFT language. In short, the basic requirement of unitarity in RFT can be formulated as a certain property of the action of the RFT Hamiltonian on the projectile and target wave functions.
Both these wave functions are constructed as superpositions of (appropriate) multidipole "Fock" states, the structure directly inherited from QCD. The coefficients of the JHEP08(2016)031 multi-dipole states, both in the projectile and target have the meaning of probabilities and hence each has to be positive and smaller than one. When acting on eitherthe projectile or the target, a unitary RFT Hamiltonian has to preserve this property. This has to hold for all projectile/target states belonging to the corresponding Hilbert spaces.
We will show that the above requirement of unitarity is not satisfied by the action of the Braun Hamiltonian on either the projectile or the target wave function. The Balitsky-Kovchegov (BK) evolution [33,34] is partly unitary, in the sense that it unitarily evolves the projectile wave function. However its action on the target wave function strongly violates unitarity. While we will not discuss this in any detail, it is clear that the same conclusions hold beyond the large N C approximation, and thus both the JIMWLK [36][37][38][39][40][41][42][43][44] and the KLWMIJ [46] Hamiltonians violate unitarity as well.
Certain problems with t-channel unitarity in the BK evolution have been already noticed a while ago in [57]. Those were believed to have been cured by inclusion of Pomeron loops along the lines of Braun's construction [30-32, 49, 58-66]. Our present analysis shows that problems with unitarity in the current RFT approaches run deeper. Although the simple prescription a la Braun is likely sufficient to restore the t-channel untarity of the BK evolution, the s-channel unitarity is violated to some degree by all currrently available implementations of high energy evolution, including the Braun version of the BFKL Pomeron calculus.
We have made an attempt to find a modified RFT Hamiltonian which implements the unitarity conditions, and also reproduces the JIMWLK and KLWMIJ evolution in appropriate limits. This attempt was so far unsuccessful in the context of the realistic 2+1 dimensional RFT. An analogous program however can be explicitly followed through in a toy model with zero transverse dimensions [67][68][69][70][71][72][73][74][75][76][77][78][79]. The standard zero dimensional toy model with triple Pomeron vertex shares the paradoxical features of the Braun theory. In the context of a zero dimensional toy model we were able to construct a modified Lagrangian, which satisfies the zero dimensional analog of the QCD unitarity conditions. This model also has the JIMWLK (or rather its BK limit [33,34]) in the appropriate kinematics.We were also able to find explicit solutions for the evolution generated by this theory, and verify that it is free from paradoxes mentioned above.
The plan of this paper is the following. In section 2 we recap the formulation of high energy evolution in the CGC approach. We also provide a path integral formulation of the calculation of the scattering amplitude, and demonstrate that in the appropriate limit it reproduces the Braun Lagrangian. We also recap the peculiarities of the high energy evolution generated by this Lagrangian.
In section 3 we confirm this strange behavior by considering small fluctuation analysis around fixed points of the Braun Lagrangian.
In section 4 we shift our attention to the zero dimensional model [67-73, 75-77, 79] in order to demonstrate explicitly that it exhibits a similar paradoxical behavior. In this context we also demonstrate in a simple and straightforward way that the evolution of the zero dimensional analog of the Braun model as well as the JIMWLK model is not unitary. Of course this statement has to be taken with a grain of salt. There is no fundamental field theory for which this model can serve as an effective high energy limit. However the JHEP08(2016)031 formal structure of the model is very similar to that of a realistic high energy QCD RFT. Thus in this context we can explore the formal analog of the QCD unitarity requirement in order to understand later its implementation in the realistic QCD RFT.
In section 5 we construct a modification of the zero dimensional toy model which satisfies the unitarity requirements. We show that this model agrees with the JIMWLK and Braun Lagrangians in appropriate limits. We provide analytic solutions for the modified unitary model, and show that this evolution is devoid of the worrisome features mentioned above.
In section 6 we return to the 2+1 dimensional QCD RFT. We show that the Braun and JIMWLK evolutions are non-unitary in this realistic context. We also discuss difficulties we face in trying to follow through the program of constructing a unitary evolution in this dimensionality.
Finally section 7 contains a short discussion of our results.

Pomeron path integral from the CGC formalism
Our goal in this section is to derive a path integral representation for the scattering amplitude starting with the expressions derived in the CGC formalism in [80]. The main motivation for this reformulation is to make direct contact with the formulation of the RFT by Braun [30][31][32].

The scattering amplitude
In the CGC formalism the scattering of a fast moving projectile on a hadronic target is given by the expression Here ρ a (x) is the color charge density of the projectile, α a T (x) is the color field of the target, and R and S are defined as with the projectile color field α a (x) determined by the projectile color charge density ρ a (x) via solution of the static Yang-Mills equations. The operator R is the "dual Wilson line". An insertion of a factor R in the amplitude eq. (2.1) is equivalent to appearance of an extra eikonal scattering factor associated with an additional parton. In this sense R creates an additional parton in the projectile wave function. The Wilson line S involves the projectile color field and has the meaning of the eikonal s-matrix of a target parton that scatters on the projectile. Here we have denoted the functional Fourier transform ofW In this paper we will adhere to the dipole model framework [7][8][9]. We therefore assume, that all the observables can be written in terms of dipoles only, in which case, neglecting the possible contribution of Odderon, the two basic elements of our calculation are the Pomeron and its dual,

JHEP08(2016)031
The integral over the charge density ρ in eq. (2.1) can be replaced by the integral over P . In principle this change of variables involves a Jacobian, but it is inessential to our discussion and we will neglect it in the following. Thus in the dipole model limit we have The structure of the weight functions W P and W T is crucially important for the subsequent discussion of unitarity. This structure has been discussed in detail [80]. The presence of a physical dipole in the projectile wave function corresponds to a factor d(x, y) ≡ 1 − P (x, y) in W P . Thus for a wave function that contains a distribution of dipole configurations (numbers and positions), the projectile weight function has the form The functions F n ({x,x}) are probability densities, and therefore are nonnegative definite F n ({x,x}) ≥ 0. Similarly, a dipole in the target wave function carries a factord(x, y) Considered as operators on the space of functionals W , the objects P andP have nontrivial commutation relations. In principle those are directly calculable from the definitions eq. (2.2), but this is not a trivial calculation. In the literature these commutation relations are usually approximated by those calculated in the dilute regime. In this regime, where any projectile dipole scatters only on a single target dipole (and vice versa), we can approximateP by [35]

JHEP08(2016)031
where γ(x, y; u, v) is the Born level scattering amplitude of a dipole (x, y) on the dipole (u, v) γ(xy, uv) = α 2 The function γ satisfies With these commutation relations eqs. (2.10), (2.11) the interpretation of the calculation of the scattering amplitude in eq. (2.4) is rather neat and intuitive. Moving one operatorP (x, y) from W T through W P one kills one of the operators P (u, v) and instead acquires a factor −γ(x, y; u, v), which is the dipole-dipole scattering amplitude. The original Pomeron P (x, y) vanishes once it arrives next to the δ(ρ). The net result is that the Pomerons P (x, y) andP (u, v) leave behind a factor −γ(x, y; u, v) and disappear from the rest of the calculation, in accordance with the approximation that any dipole of the target can only scatter on one dipole of the projectile, and after doing so does not participate in any further scatterings. This clearly corresponds to dilute limit where a given projectile dipole can meet at most one target dipole while traversing the target.
In sections 5 and 6 we will discuss the modification of the commutation relation between P andP and an associated interpretation in terms of dipole-dipole scattering.
Eq. (2.4) defines the scattering matrix at some initial rapidity. The S-matrix evolved by the rapidity Y is given by (2.14) Eq. (2.14), does not presuppose the commutation relation eq. (2.9) and remains valid if P is a more complicated function of the conjugate Pomeron P † , e.g. of the type we will discuss in subsequent sections. The Hamiltonian H RF T is known in two limits -for dilutedense [46] and dense-dilute [36][37][38][39][40][41][42][43][44] situation. Another version of H RF T was suggested by Braun [30][31][32] as appropriate to description of scattering of two nuclei. We will discuss explicitly these Hamiltonians later, but for now our goal is to derive a general path integral representation for the S-matrix eq. (2.14).

The path integral representation
First let us note that the expression eq. (2.14) is a multidimensional generalization of a "quantum mechanical" amplitude of the general form Using the exponential representation of the δ-function, and the fact that the only non-vanishing contributions come from terms where all derivatives in W 1 (p) act on this exponential, we can write

JHEP08(2016)031
With the usual trick of inserting resolution of identity at intermediate "times" (a.k.a. "rapidities), this can be written as the integral over trajectories with somewhat unusual boundary conditions Returning to the RFT, and using the correspondence which follows from the commutation relation eq. (2.10), we can write 1 (2.19) Throughout this paper, we will be focussing on several Pomeron Hamiltonians. The first one is the BK Hamiltonian, which is a dipole/large N C version of the KLWMIJ Hamiltonian: where K(x, y|z) = (x−y) 2 (x−z) 2 (y−z) 2 is the BFKL kernel in the dipole form. The dual version of this Hamiltonian, which we will refer to as the KB Hamiltonian, is a dipole/large N C version of the JIMWLK K(x, y|z) P † (x, z) +P † (z, y) −P † (x, y) −P † (x, z)P † (z, y) P (x, y) (2.21) Note that to write eqs. (2.20), (2.21) we did not have to assume a specific relation between P and P † (or P andP † ). The third interesting Hamiltonian was written by Braun in [30][31][32], and it explicitly assumes the dilute regime commutation relations eq. (2.9).
In order to write (2.19) in the form presented in [30][31][32], we introduce a "linearized" version of the Pomeron by where the last approximate equality holds for small P , or dilute projectile limit. In the same approximation

JHEP08(2016)031
We then have As for the weight functions W , it follows from our previous discussion and in particular eqs. (2.5), (2.6) that they can be expressed in term of Φ andΦ. For a projectile and a target with fixed numbers of dipoles at given coordinates we have where we have assumed that P andP are both small at initial rapidity. The "currents" J P (x, y) and J T (x, y) are simply the number density of the dipole in the projectile and target at points (x, y) respectively. These weight functions can be traded for fixed boundary conditions on the Pomeron fields. Differentiating the action with respect to Φ gives the equation of motion forΦ with the source term J P (x, y)δ(η − Y ). Similarly the equation of motion for Φ acquires the source term J T (x, y)δ(η). Integrating as usual the appropriate equation of motion across the appropriate boundary one finds that the presence of the source terms is equivalent to imposing the boundary conditions These boundary conditions are equivalent to specifying the Born amplitude for the dipole scattering on the target and the projectile prior to rapidity evolution.

Peculiarities of the Braun evolution
In this section we discuss qualitatively the nature of solutions to the equations of motion generated by the Braun Hamiltonian. The classical equations of motion are given by If the values of φ andφ are not chosen to be exactly the "fixed point" values it is obviously impossible for the solution to reach any of the fixed points at all rapidities. However one might expect that if the evolution interval in rapidity is very large, Y → ∞, the solution will be arbitrarily close to one of the fixed point values for large interval of intermediate rapidities η of the length of order Y away from the end points. This expectation may be too naive, and in fact we will see that in the zero dimensional toy model it is not strictly satisfied. However one certainly does expect that starting from the end points and moving towards the midpoint of the rapidity interval, the solution will develop towards one of the attractive fixed points, even though it may not quite reach it.
Since the interpretation of Φ andΦ is that of the scattering amplitude of an external dipole on the target and the projectile respectively, the point (0, 0) is the vacuum fixed point, where the scattering amplitude on both the target and the projectile vanishes. It is a repulsive fixed point, meaning that the solution of equations of motion departs from it with rapidity, given an initial condition which is not exactly zero.
The point (1, 1) corresponds to the dense-dense limit, where both the target and the projectile are black. Since one expects both the target and the projectile states to become dense as a result of the evolution, one expects this to be an attractive fixed point. In other words we expect that for boundary conditions φ = 0 andφ = 0 and for a very large rapidity interval Y , the solution will be approaching the point Φ(η) = 1 andΦ(η) = 1 towards the middle of the rapidity interval 0 < η < Y .
Finally we will refer to the points (0, 1) and (1, 0) as the "BK" fixed points, since they correspond to the situation where one of the colliding objects is dense and one is dilute. Naively we expect these two fixed points to be repulsive, albeit not as strongly repulsive as (0, 0). In other words, if the boundary conditions are not too far from these values, e.g. Φ(0) = ,Φ(Y ) = 1 − , for an intermediate range of Y the solution will be close to the point (0, 1) for most intermediate rapidities η. However if Y is increased to sufficiently large value, the solution for these boundary conditions will eventually flow towards (1,1) at intermediate values of Y η 0.

JHEP08(2016)031
Our goal in this section is to determine whether our intuitive expectation on the nature of the fixed points is born out by the equations of motion of the Braun model. In the following we concentrate on the points A = (1, 1) and B = (1, 0).

Point A
First, let us take φ andφ both close to unity, which means that already at initial rapidity both the projectile and the target are dense objects (nuclei). It is then reasonable to expect that Φ andΦ stay close to unity in the whole rapidity interval 0 < η < Y .
We can then write down the linear equation for the deviation of Φ andΦ from unity, or after obvious cancellations It is more transparent to multiply these equations by the operator L and write them as equations for n(xy) = L(xy)∆(xy) andn(xy) = L(xy)∆(xy). The physical meaning of n (similarlyn) is that of the logarithmic dipole density. Using the fact that the BFKL kernel commutes with the operator L we obtain: All the kernels on the r.h.s. of these equations are positive, and thus we conclude that ∂ ∂ η n(x, y; η) < 0 and ∂ ∂ ηn (x, y; η) > 0. This means that Φ approaches unity as the rapidity η increases, whileΦ approaches unity as η decreases. Thus the fixed point A is attractive, namely if we start with initial conditions where Φ is close to unity at η = 0 andΦ is close to unity at η = Y , and Y is large enough, then at all values of η the solution will be close to the fixed point. This is in accordance with our naive expectation.

Point B
Now consider point B. Assume that our initial conditions fix Φ to be close to one at η = 0, but fixΦ to be small at η = Y . Denoting 1−Φ ≡ ∆,Φ ≡∆, we have the small fluctuation equations as The first equation, as before says that Φ approaches unity at Y > 0.
The second equation when rewritten in terms ofn reads This shows thatn increases towards positive rapidities. If it starts off as small at η = Y it becomes even smaller at η < Y and thus this fixed point also is attractive. This simplified analysis suggests that both A and B are attractive fixed points, and thus depending on the initial conditions, the system flows to one or the other at intermediate rapidities. This result is surprising and goes against our intuition. It suggests that for a physical situation where the dense-dilute scattering process is evolved to large rapidities, the dilute object never gets dense.
This interpretation of the behavior we have just found is not quite adequate, as we will discuss in the next sections. Nevertheless this behavior is counter intuitive. However our analysis is incomplete, since we have been a little cavalier about the dependence of the Pomerons on transverse coordinates. Although in the strict sense the points A and B are indeed all the fixed points of the flow, one never starts with initial condition which is close to saturation at all values of dipole size. A more careful analysis should allow for existence of a finite saturation momentum. It is logically possible that accounting for finite Q s will change the attractive nature of the point B. We will now perform this analysis.

Point B -finite Q s
Let us consider again vicinity of the point B. Let us assume that at any rapidity η, . We will still assume thatΦ is small for all dipole sizes. The value of the saturation scale Q s depends on η.
Consider eq. (3.1) for the Pomeron Φ. Here the contribution of the second line in eq. (3.1) is always second order in smallness, sinceΦ is small and L annihilates the constant part of Φ. Thus the equation of motion in this approximation becomes the BK equation whose behavior is well understood.
Consider first small external dipole sizes x − y < Q −1 s (η). The contribution of the nonlinear term from large emitted dipole sizes z > Q −1 s cancels half of the contribution of JHEP08(2016)031 the real part of the BFKL kernel, while the region of small z leaves the contribution of full BFKL equation If we were to neglect the last term, the solution would be just that of the BFKL equation, namely exponentially growing with rapidity. The last term is also positive, and thus speeds up the evolution slightly. This term however is only important when the Pomeron is in the color transparency regime. Recall that the BFKL kernel decreases as z −4 at large z. Thus the contribution of the last term is proportional to α s |x − y| 2 Q 2 s = α 3 s |x − y| 2 µ 2 , where µ 2 is the gluon density. The first term is the BFKL equation with gluon emissions limited to the short distance. Its contribution can be estimated approximately as with ω a number of order unity. Eq. (3.12) then effectively reads Thus once the Pomeron surpasses its color transparency limit, the second term is negligible and the evolution is dominated by the first term which is just the BFKL evolution. For large dipoles, |x − y| > Q −1 s (η) the nonlinear term in eq. (3.1) cancels the contribution form the real part of BFKL kernel in the region |z − x| > Q −1 s (η) and |z − y| > Q −1 s (η). In this large z region we can write as before Φ(x, z) = 1−∆(x, z) etc. The only contribution to the r.h.s. of eq. (3.1) in this region is then . Substituting this into the r.h.s. of eq. (3.1) we find complete cancellation to linear order in ∆. The same happens in the second region of small z. The only leftover is the virtual term integrated over the remainder of the space |z − x| > Q −1 Since the integral of the virtual term is cut off in the infrared in the range |z − x| ∼ |x − y|, the equation for small fluctuations in the saturation region becomes This is of course, the standard Levin-Tuchin argument [81]. It shows that as the rapidity increases, the Pomeron approaches saturation. At small values it grows toward saturation JHEP08(2016)031 according to BFKL, while close to saturation, it continues to grow albeit slowly. Since Q s (η) grows with η, more and more dipole sizes are saturated as rapidity increases. The more interesting and problematic equation is the one for∆ ≡Φ (3.2). Our main interest here is to see whether allowing for a finite saturation momentum can reverse the flow ofΦ and somehow through a back door cause it to grow towards the smaller values of rapidity η < Y .
We consider the initial condition whereΦ is small at η = Y . It's saturation momentum in this rapidity range is vanishing. However in the evolution we have to account for the effect of the finite saturation momentum of Φ. Taking this into account we find that the contribution to the last two terms in eq. (3.10) is restricted to the integration region |z − x | > Q −1 s (η) and |z − η| > Q −1 s (η) respectively, since in the rest of the domain this contribution is quadratic in ∆∆.
Consider first large external dipoles |x−y| > Q −1 s . In this regime (commuting as before the operator L −1 with the BFKL kernel), the non-linear term cancels the contribution of the real BFKL kernel except in the region |z − y| < Q −1 s in the first term of eq. (3.10) and |z − x| < Q −1 s in the second term of eq. (3.10). These leftovers are almost cancelled by the appropriate part of the virtual integral. The remainder depends on the small size dipoles, so that it plays the role of a source in the equation: The only difference between this equation and eq. (3.11) is the last line. Without this term the densityn decreases towards small η, as discussed above. The source term itself is positive and thus potentially could change this behavior. The equation can be re-written in the form: To understand whether the source term can have an important effect, we have to consider the evolution of small dipoles as well.
The evolution of the small dipoles |x − y| < Q −1 s (η) is given by K (x, y|z)n (z, y, η) Here the last line appears as the source term due to coupling of large dipoles. The contribution of large dipoles with sizes greater than the inverse saturation momentum is absent in the r.h.s. of eq. (3.19), since the last two terms cancel the contribution of those JHEP08(2016)031 dipoles to the BFKL kernel. The equation forn(x − y) for small dipoles therefore is not sourced by large dipoles.

−
∂n(x, y; η) ∂η =ᾱ It is simple to solve this equation for a certain set of initial conditions. Let us consider a situation when the unevolved projectile is a single dipole of the size R 1 while the target is dense R 1 > Q −1 s (0), so that the scattering amplitude on the target is close to unity at all rapidities. For this situation the initial condition for eqs. (3.17), (3.19) is Since the initial dipole size R 1 never gets inside the saturation radius throughout the evolution, the small dipolen satisfies at all rapidities a homogeneous equation with the vanishing initial condition. Thusn(x − y; η) = 0 as long as (x − y) 2 Q 2 s < 1. This also means that there is no source term coming from small dipoles in the equation for large dipoles, and the solution forn(x − y) of arbitrary size is Although the source term can be important for other initial conditions, it is obvious from the previous discussion that at least for some physically reasonable initial conditions the introduction of finite Q s does not change the results of the previous subsection. Namely, it is indeed true that the BK fixed point (1, 0) is the attractive fixed point of the Braun Hamiltonian [56]. In principle, one could perform a similar more refined analysis of the stability of the point (1, 1) allowing for finite Q s of the projectile and the target. We will not do it here, as the paradoxical nature of the evolution generated by the Braun Hamiltonian is already clear.

What does it mean?
We have established that in the classical approximation to Braun evolution, the Pomeron Φ(η) decreases towards small values of η, if Φ(η) is close to unity. This behavior is counterintuitive. NaivelyΦ(η) has the meaning of the scattering amplitude of a dipole on the projectile wave function at rapidity Y − η, and so this seems to suggest that the projectile becomes more transparent to dipoles that have higher energy, if the target is black.
We would now like to be a little less naive and understand more formally what is the meaning of this behavior. Solving the classical equations of motion forP is the classical approximation to calculating the rapidity dependent average P (η) . Consider again the correspondence between the CGC expression eq. (2.14) and its path integral representation eq. (2.24). It is clear from this correspondence that the rapidity dependent average ofP JHEP08(2016)031 in the CGC formulation is given by the expression where W η T [P ] is the target wave function evolved through the rapidity interval η. This is the scattering matrix of the projectile on an object which is obtained by evolving the target by rapidity η, adding to it one extra dipole, and then evolving the resulting system by rapidity Y − η. What does one expect the η dependence of such a scattering matrix to be? Clearly, adding an extra dipole towards the end of the evolution of the target should be less efficient in making the target black than adding it earlier in the evolution. If one adds a dipole early on, it should contribute to subsequent evolution and lead eventually to relatively more dipoles in the wave function, since the QCD evolution always increases the number of physical dipoles. Thus we expect that 1 −P (η) should increase with η, and thereforeP (η) should decrease with η, or equivalentlyP (η) should increase to smaller values of η.
As we saw earlier, the behavior ofP in the Braun evolution is opposite. There are two possible interpretations of such behavior. One is that the target wave function becomes less dense with evolution. In this case the additional dipole is "bleached" by further evolution and contributes little to the scattering amplitude. Physically such behavior is of course completely unacceptable. The other possibility is that the target does become denser, but the evolution is so violent that dipoles disappear from its wave function by physically merging with each other. This possibility may be more palatable a priori, however it also is not a part of QCD dynamics. In the leading logarithmic approximation the QCD evolution produces new gluons, and therefore dipoles and never annihilates partons that already exist in the wave function. QCD saturation is the statement that the rate of this growth decreases with the density of the target, but it never completely vanishes and certainly does not become negative. Thus the behavior of the solutions of the Braun equations indeed violates QCD expectations.
One could wonder if perhaps this means that the classical approximation to the path integral, which yields this behavior is violated in the dense-dilute regime. It is in principle possible, since the applicability of classical approximation is determined by the magnitude of the sources J T and J P , and in the dense-dilute limit one of the sources J P is small. In this case one should be able to see that the loop corrections to the classical approximation are large. However it has been shown in [82,83] that the presence of large Φ leads to appearance of a kind of "mass" in the Pomeron propagator, and in fact suppresses the loop corrections. The culprit therefore seems to be the Braun Hamiltonian itself, rather that the classical approximation to the evolution.
We note that this behavior is not unique to the Braun Hamiltonian. In particular we can repeat the same analysis in the framework of the Balitsky-Kovchegov equation, or the dipole model approximation to the JIMWLK Hamiltonian. The result is exactly the same. The evolution of a dense target seems to "bleach" the scattering amplitude of an extra JHEP08(2016)031 dipole that is added to it. The BK (and JIMWLK) evolution leads to the fixed point (1, 0) at large rapidities.
This strange behavior leads one to suspect that not all is right with unitarity in the Braun and BK evolution. In particular the evolution of the dense state in any of these frameworks may be violating the QCD unitarity. In the next section we will consider this question in the context of a toy model with no transverse dimensions. We will return to the realistic case later.

Playing with toys I: trouble in the toy world
In this section we consider a set of toy models in the framework of the zero transverse dimensional reduction of the RFT [67][68][69][70][71][72][73][74][75][76][77][78][79]. Such models have long been used as a simplified setup for qualitatively understanding of the high energy behavior of QCD. Like in the previous sections we first formulate these models in the framework of the zero dimensional analog of the CGC formalism, and explore their properties. The scattering matrix of the projectile consisting of m dipoles on the target consisting ofn dipoles is given in analogy with the real QCD case by This amplitude is evolved in energy according to In the following we will consider several toy Hamiltonians.

The BK evolution
The zero dimensional analog of the BK evolution [33,34] is given by the Hamiltonian 2 As before we take P andP to have the dilute limit algebra, such that The constant γ is the zero dimensional proxy for the dipole-dipole scattering probability. The scattering matrix can be calculated explicitly (we assume m + 1 <n) In particular 1|n = 1 −nγ (4.6)

JHEP08(2016)031
It is clear from eq. (4.6) that our current formulation is restricted to number of dipoles n < 1/γ, since otherwise the single dipole S-matrix becomes negative. The reason for this unphysical behavior is the dilute limit commutation relation we adopted for P and P . As we will see later, with correct commutation relation the problem does not arise. Nevertheless as long asn < 1/γ we can continue the present analysis.
Our aim now is to compare P (η) for two different values of rapidity. For simplicity we will take the rapidity interval to be infinitesimally small, and will simply insertP into the matrix element either before or after this short evolution interval. Also for simplicity we take the projectile to contain a single dipole, although the calculation can be easily generalized. Thus we are interested in 1 −P P ≡ 1|(1 −P )e H∆ |n ; 1 −P T ≡ 1|e H∆ (1 −P )|n (4.7) As discussed in the previous section, the two quantities have simple physical meaning. The first one corresponds to the evolution of the target wave function with subsequent insertion of an additional dipole, prior to scattering on the projectile. In the second one we insert an extra dipole into the target wave function, then evolve the combined target+dipole, and then scatter it on the projectile. The physical expectation is that inserting the extra dipole closer to the target will produce a blacker target. Thus we expect 1 −P P > 1 −P T ? (4.8) As we will see, this expectation is fulfilled whenn ∼ 1. However the evolution produces the opposite result whennγ ∼ 1, that is for dense target.

JHEP08(2016)031
Thus for smalln eq. (4.8) is satisfied, however forn > 1 2γ the situation is reversed. The toy model thus behaves in the similar way to the model with two transverse dimensions.
To understand the origin of this behavior consider the infinitesimal evolution of the projectile and target wave functions with the Hamiltonian H: e ∆H |n = (1 + ∆n)|n − ∆n[1 + γ(n − 1)]|n − 1 + ∆γn(n − 1)|n − 2 (4.17) Note the sea of difference between the two expressions. Recall that the coefficient in front of an n-dipole state has the meaning of probability to find this number of dipoles in the the wave function. The projectile evolution is unitary: all the probabilities in the evolved state remain positive and smaller than unity, and the sum of the probabilities adds up to unity. On the other hand the target evolution is non-unitary. The probability to find the initial state |n after a short interval of evolution exceeds unity, while the probability to find a state with one less particle is negative. The coefficients still sum to unity like for the projectile, but clearly the target evolution violates unitarity.
We stress that the probabilities in question are probabilitites to find physical dipoles in the wave function of the evolved hadronic state. Thus the negativity of probabilities violates the s -channnel unitarity. This violation is not directly seen in the calculation of the diagonal matrix element of the S-matrix, since by construction this matrix element is a real number smaller than one. However it is clear that if we were to consider more exclusive observabes, the negativity of probabilities would show up as unphysical values for some observable. Although a detailed study of this question is outside the scope of this paper, it is easy to give an example of such an observable. Consider a target containing n dipoles, which scatters on a dense projectile. If the projectile is very dense, all the dipoles in the target wave function will be scattered into the final state. Thus the probability for producing n dipoles in the final state will be equal to the probability of finding n dipoles in the wave function. At a slightly higher energy than the initial one, the probability of producing n − 1 dipoles in the final state will be negative. Strictly speaking those are of course "toy dipoles" in the "toy hadron", but the essence of the argument is the same in real QCD.
It is interesting to examine more carefully the evolved target side wave functions that enter in the calculations in eq. (4.14) at the rapidity they scatter on the projectile dipole. Calculating the average number of dipoles in the two wave function we find N T =n + 1 + ∆(n + 1) − ∆γ(n + 1)n (4.20) N P =n + 1 + ∆n − ∆γ(n − 1)n (4.21)

JHEP08(2016)031
This indeed displays the poignant feature discussed above, namely at smalln the average number of dipoles in the wave function is larger if an extra dipole is added before evolution, while at largen > 1/2γ the situation is reversed. The reason for negative difference in eq. (4.20) is obvious. It appears because the negative probability of the |n−1 contribution in eq. (4.17) grows withn faster than the positive probability of the |n contribution.
Thus we see that the nonunitarity of the BK evolution of the target wave function is indeed the reason for the counter intuitive behavior ofP with rapidity.
To summarize, we have shown that although the (zero dimensional) BK-JIMWLK evolution of the projectile wave function preserves QCD unitarity, the same evolution when viewed as evolution of the target wave function is non-unitary.

The Braun Hamiltonian
Next consider the analog of the Braun Hamiltonian We pose the same question: does this Hamiltonian generate a unitary evolution? To answer this we consider, as before e ∆H B |n ≈ (1 + ∆H)|n = (1 − ∆n)|n + ∆n|n + 1 − ∆γn(n − 1)]|n − 1 + ∆γn(n − 1)|n − 2 (4.24) This is somewhat more satisfactory than eq. (4.17), since the violation of unitarity is O(γ) and is small for smalln . However the coefficient of the term |n is still negative, and becomes large parametrically long before the saturation limit is reached. Alarmingly, since the Braun evolution is symmetric between the target and the projectile, the projectile evolution now is also non-unitary and involves negative probabilities.
We note, that the Braun eq. (4.23) has been considered in the past from the point of view of the reaction-diffusion process (RDP) [84]. Ref. [84] indeed made it explicit that this evolution corresponds to a non-unitary RDP that involves negative emission probabilities. The RDP emission probabilities are however distinct from the QCD probabilities and in fact not related to them in a simple obvious way. Thus the violation of unitary we discuss here is distinct from, and not obviously related to the nonunitarity of the appropriate RDP.
There may be more than one problem in the previous models. In particular we have seen that the commutator we have postulated between P andP can only be used for a target with small enough number of dipoles, otherwise even without any evolution the S-matrix is non-unitary. In particular 1|n < 0 for large enoughn > 1/γ. One could perhaps wonder if this deficiency is to blame for the nonunitarity of the evolution as well.
In the rest of this this section we will rectify this deficiency and show how to define the correct commutation relation. We will also show that even with the redefined commutation relation, the BK and Braun Hamiltonians lead to non-unitary evolution.

Are the commutators to blame?
Our postulated commutation relation does not allow for multiple scattering corrections when a single dipole of the projectile scatters on several dipoles of the target. Clearly the JHEP08(2016)031 correct formula for scattering of one dipole onn dipoles should be 1|n =n since this expression correctly accounts for multiple scattering corrections. The algebra of P andP should be such that this result follows from the definition of the amplitude eq. (4.1).
A simple way to achieve this is to modify the relation between P andP as follows This is better, but still not good enough. In particular it does not allow two dipoles of the projectile to scatter on the same dipole of the target, since the first factor 1 − P by differentiation simply kills the particular target dipole, and subsequent scatterings on it are not possible. The propagation of the projectile dipole should be "non demolition", in the sense that after moving the factor 1 − P throughP , the factorP should not disappear from the wave function W T . Therefore a more reasonable representation is However eq. (4.27) is not quite adequate either. According to it the propagation of the projectile dipole does not destroy any target dipoles, but the projectile dipole itself disappears after propagation, and this is not right. None of the above problems arise if the P andP algebra is taken to be the following This ensures that moving one projectile dipole throughn target dipoles give the correct factor (1 − γ)n that includes all multiple scattering corrections, while all the dipoles remain intact, and can subsequently scatter on additional projectile or target dipoles. For small γ and in the regime where P andP are small themselves, we obtain [P,P ] = −γ + . . . (4.29) consistently with our original expression (2.10), (2.11). Note that the algebra eq. (4.28) is equivalent to the following representation In the calculation of an amplitude of the type of eq. (4.1), once all the factors of 1 −P are commuted through to the left, in any matrix elementP hits the δ-function and thus vanishes. The remaining factors of (1−P ) also turn to unity, since a factor of Φ is equivalent to a derivative acting on the δ-function, and when integrated overP vanishes.
With the new algebra we have which is a simple and intuitive result: the s-matrix of dipole-dipole scattering to the power of the number of dipole pairs that scatter.

JHEP08(2016)031
We stress that the modification of the Pomeron algebra is not a matter of choice, but is necessary to obtain the amplitude eq. (4.31), which is unitary for arbitrary numbers of colliding dipoles. However the question of the unitarity of the evolution is a completely separate one. We will now reexamine the BK and Braun evolutions with the modified Pomeron algebra.

BK evolution revisited: the Hamiltonian with modified commutators
We start with the BK Hamiltonian defined in eq. (4.3). As before we ask if evolution by the infinitesimal rapidity interval preserves the probabilistic interpretation of the initial wave function.
This result is surprising: the evolution of the target is now unitary, but of the projectile is not. The evolution on the target wave function looks reasonable. Whenn is small, it is identical with the BFKL evolution. For largen it exhibits very strong saturation due to suppression with the factor (1 −γ)n, so that at largen the evolution is super slow. This is a little disturbing, but does not seem fatal. However the projectile now evolves nonunitarity, and thus we expect the same type of trouble as found in the previous subsection.
Let us see how this reflects on the behavior of P andP . Calculating P we find This difference is always negative, consistently with our logic, even though the evolution as we saw is nonunitary. Now for P : This has the same behavior as before. For smalln this difference is positive, while for n > − ln(2−γ) ln(1−γ) it changes sign, and thus it again manifests nonunitarity of the evolution.

The Braun Hamiltonian with modified commutators
Next let us examine the evolution generated by the Braun Hamiltonian. The analog of the original Braun Hamiltonian eq. (2.27) is with Φ = ln(1 − P );Φ = ln(1 −P ). The variables Φ andΦ are canonically conjugate. The modified commutation relation between P andP however enters in the calculation of matrix elements, as the projectile and target wave functios carry factors of (1 − P ) n ; (1 −P )n, see eqs. (2.5), (2.6). The action of this Hamiltonian is obviously nonunitary, since a basic necessary condition is that the coefficients in the expansion of H in powers of (1 − P ) be finite. The logarithmic factors in eq. (4.36) obviously yield infinite expasion coefficients. However using instead the form eq. (4.23), which is self dual and reduces to eq. (4.36) in the dilute limit solves this problem. This form of the Braun Hamiltonian has a chance to be unitary, and we will concentrate on this question. We rewrite the Braun Hamiltonian eq. (4.23) in a more convenient form: The action on the projectile and the target is obviously symmetric, as the Hamiltonian is self dual under the transformation P →P .
This is a nasty surprise. Now the evolution of both, projectile and target is non-unitary.
In fact the lack of unitarity is there for arbitrary number of dipoles m andn. We may hope that modifying the Braun Hamiltonian with an extraP 2 P 2 term could improve the situation. However, it does not bring about complete redemption. Consider (4.40) Now we obtain: Unfortunately this is as non-unitary as the BK evolution with the slight modification of the approach to saturation. It is not difficult to show that the nonunitarity cannot be cured by adding the four Pomeron vertex with any coefficient.

JHEP08(2016)031
5 Playing with toys II: making the toy world a better place It may seem that our modification of the commutation relations was in vain as it did not solve the problem of non-unitary evolution. However, as it happens often, good deeds get rewarded. In this section using the correct commutation relations we will be able to find a Hamiltonian which has the correct dense-dilute limit, is self dual and produces unitary evolution of both, the projectile and the target.

Unitarity regained
The discussion of the previous section does not necessarily mean that we are doomed to live with non-unitary evolution. There is one thing that we have so far implicitly accepted, namely the form of the BK/Braun hamiltonian in terms of the Pomeron operators. This is so even though we have not derived it directly. What is derivable from QCD is the Hamiltonian in terms of P and P † , rather than P andP . Since P † andP are only proportional to each other in the dilute limit, our use of the BK and Braun Hamiltonians away from this limit is not justifiable. We do know however, that the correct unitary Hamiltonian (if it exists) has to reduce to the H BK in the limit of smallP . We will now attempt to modify the BK Hamiltonian in a way that makes it unitary, but still reduces to the original H BK when expanded to linear order inP . First, we express the P † in terms ofP . To do this recall that P † should annihilate a dipole when acting on the wave function. Using eq. (4.30) we can write Thus the BK Hamiltonian expressed in terms of P andP is We can also conveniently write its dual (the mean field approximation to the KLWMIJ Hamiltonian) as In the above equations for simplicity we have used ln(1 − γ) ≈ −γ, since γ ∼ α 2 s 1. To define the Braun Hamiltonian we have to add these two and subtract the BFKL limit. The simplest analog of the BFKL Hamiltonian is the leading order expansion of either one of eq. (5.2) or eq. (5.3) in Φ and d/dΦ We thus can write an analog of Braun Hamiltonian as

JHEP08(2016)031
This Hamiltonian is clearly non-unitary. One does not need to perform any calculation to understand this. Our unitarity test of the projectile evolution amounts to the following simple three step procedure: 1. Act with the Hamiltonian on a monomial (1 − P ) n ; 2. Expand the result in powers of (1 − P ); 3. Check that the coefficients of all terms (1 − P ) m ; m = n are positive, and the coefficient of (1 − P ) n is negative.
For the target the same procedure is applied to (1 −P ) n . This set of conditions can be formulated as the following requirements on the Hamiltonian. Write the Hamiltonian as a function of d = 1 − P andd = 1 −P ; H(d,d). When the hamiltonian is commuted throughd n to the right, each operator d turns into one-on-n scattering amplitude, a positive number smaller than one, which we can also denote as d. Thus our requirements can be written as It is obvious that neither H 1 B , nor H BK nor H KB is unitary, since they all contain logarithmic factors. Thus step 2 in our procedure fails, as it gives infinite coefficients. Equivalently, the derivatives of the Hamiltonian atd = 0 diverge.
However the proposal eq. (5.5) for the Braun Hamiltonian is not unique. It was written on the basis of two requirements: it should be self dual under P ↔P ; and for smallP (P ) it should reduce to H BK (H KB ). It is in fact possible to write down a Hamiltonian that satisfies these requirements, as well as the unitarity constraint: where U T M stands for "Unitarized Toy Model". The fact that it is self-dual is evident. ExpandingP to linear order in P † , using eq. (5.1) leads to the BK Hamiltonian, eq. (5.2).
To check the unitarity we consider: This evolution is clearly unitary. Due to self duality, it is clear that the evolution of the projectile wave function is unitary as well. Interestingly it also exhibits the saturation behavior very similar to the one that is expected from the real QCD evolution, namely at largen, the change in the wave function is independent of the number of dipolesn. In this respect it contrasts strongly with eqs. (4.16) and (4.32), which were also unitary. In the standard BK evolution of the projectile eq. (4.16) the wave function never saturates, meaning the rate of growth of number of dipoles is proportional to the number m of dipoles JHEP08(2016)031 in the state even for very large m. This is of course the well known property of the BK evolution, where the projectile state evolves according to the perturbative dipole model and saturation of the scattering amplitude is only due to the multiple scattering effects. Eq. (4.32) on the other hand is completely different. Its evolution is "oversaturated", in the sense that for largen, the wave function does not evolve at all. Clearly this cannot be a reflection of a QCD-like dynamics. An interesting and very appealing property of the UTM Hamiltonian, is that one can arrive at it either from BK by expanding P † to leading order inP , or from KB by expandinḡ P † to leading order in P , or indeed from BFKL by using both expansions.
Having found a unitary evolution it is interesting to explore its properties. In the next subsection we provide a solution of the classical equations that follow from H U T M .

Equations of motion and the scattering amplitude
The general form of equation of motion follows from dP dη = H, P ; and dP dη = H,P (5.9) With the Hamiltonian H B we get Interestingly, although it is not obvious from the form of the Hamiltonian eq. (5.7), the evolution has the same fixed points as in two transverse dimensions (0, 0), (1, 0), (0, 1), (1, 1). Since the Hamiltonian is conserved, we havē The general solution to eq. (5.10) takes the form: where the parameters β and α should be found from the boundary conditions: One can see that for p 0 >p 0 and e (1−α)Y 1, eq. (5.13) leads to For a symmetric boundary condition p 0 =p 0 eq. (5.13) give P (0) =P (Y ) and the solution takes the form
In figure 1 we have plotted numerical solutions to eq. (5.10) that correspond to different initial conditions. The BFKL-like regime is clearly seen on figure 1-a. All the solutions clearly show that P grows towards positive rapidities, whileP grows towards negative rapidities. This is of course a direct consequence of the conservation ofP P , and thus the unitarized evolution indeed cures the peculiarity of the evolution ofP .
However we learn from these solutions that our initial expectation that for large Y at intermediate rapidities the solution should be dominated by the fixed point (1, 1) is not warranted. Although both P andP grow towards midrapidity, there is no value of rapidity at which they are simultaneously close to unity, unless it is forced by the initial conditions.
In fact, once we account for the conservation ofP P , we get a very different view of the fixed point structure of the evolution. Plugging the relation eq. (5.11) in eq. (5.10) we obtain the following equations: The fixed points (0, 0), (0, 1) and (1, 0) are not present in these equations, which means that neither one of them is reachable at α = 0. The point (1, 1) is also unreachable by the evolution for α = 1. Eq. (5.17) has only two interesting fixed points: (1, α) and (α, 1). Since for any physical initial condition P (0) > α, the asymptotics at η → ∞ is always dominated by the fixed point (P = 1,P = α), while for η → 0 the point (P = α,P = 1) is approached. Whether either one of these points is reached during the evolution to finite Y depends on the initial conditions. As illustrated on figure 1-a,b for symmetric initial condition the solution approaches very close to the fixed points at both ends, while for asymmetric initial conditions this is not the case, and only the vicinity of one fixed point is reached. In figure 1-c-f this is the point (1, α) at η → Y . There are of course mirror solutions where instead the point (α, 1) is approached at η → 0.
Another noteworthy property of the solution is, that for strongly asymmetric initial conditions, the smaller of the two Pomerons remains small essentially over the whole evolution. This is clearly seen in figure 1-d. The physical reason for that is the saturation effects in the wave function. As we have seen in eq. (5.8), when the target wave function contains many dipoles (P is close to unity), the rate of increase of the dipole number is constant and independent on the number of dipoles . Consider the dependence ofP on η.
As we have discussed in detail in the previous sections, the classical solution for 1 −P (η) is the scattering amplitude on the projectile of the target with an extra dipole inserted at rapidity η. Inserting an extra dipole at rapidities η into the target wave function in principle affects the evolution of the target wave function between rapidity η and rapidity Y , at which the target scatters on the projectile. However, since in the dense regime the JHEP08(2016)031 rate of the evolution does not depend on the number of dipoles, there is in fact almost no dependence on η as long as at that η the target is dense (P is close to unity). ThusP is a nontrivial function of η only in the rapidity interval in which P significantly differs form unity. This is clearly illustrated on figure 1-d.
Clearly the classical solutions of eq. (5.10) determine the scattering amplitude in the semiclassical regime. To see this explicitly we employ the path integral representation for the S-matrix. The UTM Pomeron Lagrangian is 18) The S-matrix is then given by In the classical approximation 3 Note that the solution forP is irrelevant for the BK amplitude, which is determined entirely by P (Y ). On the other hand the scattering amplitude in UTM does depend on P . Nevertheless the two models should be close in the regime where the BK evolution applies. An example of numerically computed amplitudes for BK and UTM with generic initial conditions is presented on figure 2 for one projectile dipole (m = 1), A = 1 − S 1n . The difference between the two amplitudes is indeed quite small, reaching the maximum of about ∼ 3.5% in the pre-saturation region. However, close to the saturation, the difference between the S-matrices of BK and UTM are quite significant, since the S-matrix itself is close to zero. In other words, unitarization significantly modifies the way the S-matrix approaches zero, i.e. the zero dimensional "Levin-Tuchin law" [81].
To illustrate this fact we plot the relative difference of the S-matrix between the two calculations in figure 3. We expect that the effect of the unitarization should be even more pronounced in less inclusive observables like particle multiplicity.   6 Two transverse dimensions?
The next natural step is to try and generalize the considerations of the previous section to the real world of two transverse dimensions. The logic of changing the commutation relations between P andP is the same as in the toy model. Thus the correct commutation relations between the Pomeron and its dual are In principle, the commutation relations (6.1) should be derivable in the large N C limit from the definitions (2.3). In the following we will assume that P (xy) = P (yx), and neglect the possible contribution of Odderon to dipole scattering. This allows the following representation P (x, y) = 1 − exp{−Φ(x, y)};P (x, y) = 1 − exp u,v γ(xy; uv) δ δΦ(u, v) (6.2)

JHEP08(2016)031
contained in the wave function. If this is the case, then the factor 1/d(x,ȳ) cancels, and the result is expandable in Taylor series. However the mutual non-locality of P † andP makes this very difficult to achieve. Perhaps some additional physical insight is needed to resolve this question. We are therefore forced to postpone this problem until better times.

Conclusions
In this paper we have examined the question of the s-channel unitarity of the QCD Reggeon Field Theory. We have shown, starting from the QCD definition of scattering amplitudes, how the requirement of unitarity of the QCD evolution should be reflected in the action of RFT Hamiltonian on the projectile and target wave functions. Our finding is that the action of the BK (and JIMWLK) Hamiltonian is unitary on the projectile wave function, but violates unitarity of the target wave function. The Braun Hamiltonian, which is a self dual extension of the BK Hamiltonian, turns out to violate unitarity of both, the projectile and the target. This unitarity violation is small in the regime where the respective Hamiltonians are applicable. However the peculiar behavior of the solutions to Braun Hamiltonian at large rapidities is closely related to this violation of unitarity. Recall that classical solutions to the Braun equations of motion bifurcate beyond some critical rapidity Y c [56]. Starting at Y c , the dependence of eitherP or P on η becomes unphysical, indicating large unitarity violation. Therefore starting from this total rapidity the Braun evolution is not trustworthy as unitarity violating effects are large.
To elucidate the unitarity considerations we have discussed toy models of RFT in zero transverse dimensions. We have found that with correct commutation relations between the Pomeron and its dual, P andP it is possible to modify the Braun evolution to make it unitary. This unitarized toy model (UTM) has many desirable properties. It is self dual, just like the Braun Hamiltonian is, reduces to BK evolution in the limit of dilute projectile and evolves both the projectile and the target wave functions in a unitary way. It also exhibits approach to saturation similar to that we expect in QCD, namely for large dipole density the rate of growth of the dipole number becomes independent of the dipole number itself.
We found analytic solution of the classical equations of motion of UTM and compared the scattering amplitude calculated in classical approximation to that of the BK model. As expected, the evolution in UTM is somewhat slower, as it takes into account the saturation effects in the projectile wave function. For dilute projectile the difference between the scattering amplitude calculated in BK and UTM models is quite small, but the pre-asymptotic behaviour differs significantly in the saturation regime.
In the two dimensional case we have provided the corrected commutation relations between P andP valid beyond the dilute limit. Unfortunately so far we were unable to find a modified unitarized version of the Braun Hamiltonian in QCD. This is left for future work.
Finally we stress that our focus in this paper was on the s-channel unitarity. For the RFT to be fully consistent, it has to be t-channel unitary as well. The BFKL and

JHEP08(2016)031
Braun Hamiltonians satisfy the t-channel unitarity constraints by construction, however their applicability is limited to scattering of dilute systems. On the other hand the BK and JIMWLK Hamiltonians are applicable to scattering processes involving one dilute and one dense system, but they lack the t-channel unitarity as discussed in [57]. It appears that self duality of a Hamiltonian ensures the t-channel unitarity of the RFT. Thus in order to generalize the RFT approach to nucleus -nucleus collisions it is imperative to find a s-channel unitary and self dual extension of the BK and Braun Hamiltonians.