The Color Glass Condensate density matrix: Lindblad evolution, entanglement entropy and Wigner functional

We introduce the notion of the Color Glass Condensate (CGC) density matrix $\hat\rho$. This generalizes the concept of probability density for the distribution of the color charges in the hadronic wave function and is consistent with understanding the CGC as an effective theory after integration of part of the hadronic degrees of freedom. We derive the evolution equations for the density matrix and show that the JIMWLK evolution equation arises here as the evolution of diagonal matrix elements of $\hat\rho$ in the color charge density basis. We analyze the behavior of this density matrix under high energy evolution and show that its purity decreases with energy. We show that the evolution equation for the density matrix has the celebrated Kossakowsky-Lindblad form describing the non-unitary evolution of the density matrix of an open system. Additionally, we consider the dilute limit and demonstrate that, at large rapidity, the entanglement entropy of the density matrix grows linearly with rapidity according to $d S_e / d y = \gamma$, where $\gamma$ is the leading BFKL eigenvalue. We also discuss the evolution of $\hat\rho$ in the saturated regime and relate it to the Levin-Tuchin law and find that the entropy again grows linearly with rapidity, but at a slower rate. By analyzing the dense and dilute regimes of the full density matrix we are able to establish a duality between the regimes. Finally we introduce the Wigner functional derived from this density matrix and discuss how it can be used to determine the distribution of color currents, which may be instrumental in understanding dynamical features of QCD at high energy.


I. INTRODUCTION
High-energy hadronic collisions at RHIC and the LHC have demonstrated an unexpected collective behavior in particle production. In particular, multiple observations of the structure of final states in p-p and p-Pb collisions at the LHC, see e.g. [1][2][3][4][5], indicate a very nontrivial dynamics that leads to a correlated structure between produced particles.
The current understanding of the origin of these correlations is based on two concurrent pictures: the dominance of the collisions -final state interactions, which in today's most popular incarnation is described in terms of transport and hydrodynamics, and the dominance of the correlations in the wave functions at the initial state and very early stage of the collision. The latter option is most commonly analyzed in the framework of the Color Glass Condensate (CGC), for a review see e.g. Ref. [6][7][8][9].
The CGC approach as applied to date has one unappealing feature. It can accommodate correlations between partons either in coordinate space, or in momentum space, but it does not give one a handle to explore the correlations between the two. On the other hand, one can expect such correlations on general grounds.
The present work is motivated by the idea that one should be able to get some sense about the correlation between the profile of the charge density and the current density in the hadronic wave function at high energy. One may expect, for example, that if the scattering catches a configuration in the projectile wave function with large density or large density gradient, this configuration will naturally also have a large current, since it does not want to stay static for a long time. The large current then may translate into relatively high momenta of the particles in the wave function. That way we may be able to relate on the level of the initial hadronic state the density (number) fluctuation in the configuration space and momentum distribution of particles.
To achieve this goal we have to expand the usual CGC vocabulary and introduce a novel concept in this field of studies -the CGC density matrix. In the standard CGC approach one is only interested in the diagonal matrix elements of the density matrixρ (usually denoted by W ), which is sufficient for the calculation of a large number of observables. However, in order to study the properties we alluded to before, the off-diagonal elements of the density matrix are required. We note that a very similar problem arose in two and many gluon production at different rapidities [10] (see also [11]); in part, our current work will rely on intuition gained and guidance obtained from Ref. [10].
The plan of this paper is the following. In Sec. II we define the concept of the CGC density matrixρ. We stress that this is a completely different object than the one used to calculate the entropy of soft gluons for example in [12,13]. In Sec. III we derive the evolution of this object with rapidity. The evolution equation turns out to be a natural generalization of the JIMWLK equation, and is of the Kossakowski-Lindblad form as could have been expected on general grounds. In Sec. IV we consider the evolution ofρ in the weak and strong field regime in turn. We take a reasonable Gaussian ansatz forρ and calculate the evolution of the parameters. We also calculate the entanglement entropy associated withρ and show that is grows linearly with rapidity. The rate of growth is given by the leading BFKL exponent in the weak field case, and half of that in the saturation regime. In Sec. V we define the Wigner functional associated withρ. Finally we close in Sec. VI with a short discussion.

II. THE CGC DENSITY MATRIX
Recall that the CGC is an "effective field theory" of high energy scattering. This is a somewhat loose term, but it does in fact have a fairly precise meaning given our standard CGC calculations. The standard practice in deriving the CGC wave function and the corresponding evolution equation is to integrate out all degrees of freedom in the hadronic wave function except for the integrated color charge density j a (x ⊥ ) = j a (x ⊥ , x − )dx − . (Note that here and below we use j = j a t a to denote the color charge densities instead of commonly used ρ. This change of notation is implemented to differentiate the color charge density and the density matrix.) Thus CGC is in fact the effective field theory on the Hilbert space spanned by j. There is one subtlety here: in general, the components j a are not independent degrees of freedom, as they do not commute with each other. In the dense limit however the commutator between j's can be neglected and they can be treated as independent. We will follow this approach in the present paper and will treat j a as commuting.
Put in this way, our ambient CGC calculations are equivalent to simply integrating out a subset of quantum degrees of freedom in the hadronic wave function. Then it is natural to ask what is the reduced density matrix on the subspace of the Hilbert space spanned by the remaining quantum degrees of freedom. In particular one can ask what is the entanglement entropy of this reduced density matrix, and how it evolves with energy (for indirectly related calculations of the entanglement entropy in high energy collisions see Refs. [14][15][16][17][18][19][20][21][22]). These are the questions we will be gearing to ask in the present work.
One may wonder why it is that in all the CGC calculations to date [56] there was no apparent need to define the full density matrix. The answer is that the knowledge of the full density matrix is not always necessary. In particular, if one considers the calculation of observables which are functions of j only, and not of their canonical conjugates, it is sufficient to know the diagonal matrix elements ofρ in the basis of charge density j. Most of the observables considered so far, like observables involving only the softest gluons in the CGC wave function, are of this type.
To elaborate on this further, suppose the density matrix of the valence gluons in this particular basis is known and can be written as (we suppress the color and coordinate indexes for simplicity) We thus conclude that indeed for the operators of this type we only need to know the diagonal matrix elements of ρ. These diagonal matrix elements are encoded in the probability density functional W [j] routinely used in the CGC approach. However, this is not the only set of operators that are of interest in high energy scattering. One stark example of an operator of a different kind is the S matrix for dense-dense scattering. The eikonalŜ-matrix for scattering on a strong color field acts on the color charge density operators by rotating them by the eikonal phasê is an arbitrary element of the SU(N ) group and may be arbitrarily far away from the unit matrix. It is therefore obvious that theŜ-matrix is not a diagonal operator in the j basis, and so non-diagonal matrix elements of the density matrix must be important in its evaluation. Another example of this type of observable is the multi-gluon production probability where gluons are produced at different rapidities. Indeed, when calculating multiple gluon production where the rapidities of gluons were significantly different [10], see also Ref. [11], from the target perspective it was necessary to introduce novel weight functionals that depended on two different j's (or related to them by eikonal factors S, see below).
It is obvious that the knowledge of W [j] is not sufficient to determine the complete density matrix. In particular, in the MV model [23,24], W [j] = exp(− 4j 2 µ 2 ) (in this paper, in order to simplify equations, we we deviate from the conventional normalization of µ 2 ); this weight functional could correspond to a variety of different density matrices with very contrasting properties. One example would be This density matrix has a factorized form, and evidently corresponds to a pure state on the reduced Hilbert space. Obviously, this is not the only possibility. A priori any density matrix of the form with real functions µ, λ, A of variables x ⊥ and y ⊥ is an allowed density matrix inasmuch as it reduces to the MV model for diagonal elements and is Hermitian. There is only one restriction. The parameters µ, λ, A have to satisfy additional constraints in order forρ to have probabilistic interpretation, i.e. all eigenvalues have to be positive (the overall normalization can always be adjusted). This is equivalent to the requirement that for any positive integer n Here the determinants are in the transverse position space [57]. We also will often use the inverse functions, e.g. µ −2 (x ⊥ , y ⊥ ), defined as Given thatρ contains more information than W [j], the natural first question is how does it evolve to high energy? To start answering this question, we first point out that previously, in Ref. [25], the evolution of the density functional for two gluon production at two significantly different rapidities was derived. The problem at hand is very similar to that discussed in Ref. [25]. Although at the time of writing of Ref. [25] the newly introduced weight functional was not interpreted as a density matrix, we will show below that the evolution derived in Ref. [25] can indeed be mapped on to the evolution of the density matrix. The evolution in question is where Q is defined by where tr c denotes the trace over color indexes. Here, as usual, S is the eikonal phase matrix for scattering of a probe gluon on the wavefunction [58]. The matrix S is determined by the color charge density via In the dilute limit (small color charge density) we explicitly have (3.6) In the following we allow ourselves to denote the argument of the density matrix intermittently by either j, α or S, as all these objects are algebraically related to each other. Recall that the standard JIMWLK Hamiltonian is given in terms of Q's as Thus, Eq. (3.1) generalizes JIMWLK evolution equation [26][27][28][29][30][31][32] to the full density matrix.

A. Derivation of high energy evolution
In order to derive this equation, we will follow the same main steps as in Ref. [10]. Our discussion will be in the framework of scattering of some dilute projectile on a dense target.
Consider an observableÔ which depends only on the projectile degrees of freedom. The projectile scatters on the target, and the observable is measured in the asymptotic state long time after the scattering has taken place. The total rapidity interval in the collision is Y . We assume that the target has been boosted to rapidity Y 0 and that the operatorÔ depends on degrees of freedom between the rapidity Y 0 and Y . In other words, the operatorÔ itself does not depend on the target degrees of freedom. It is an operator in the projectile Hilbert space and as such defines an observable measured in the direction of the projectile.
Let us define the following object where |P Y −Y0 is the wave function of the projectile. This object is related to the actual observable in the scattering process via whereρ Y0 is the target density matrix we are interested in. Note that although the operatorÔ itself does not depend on the target degrees of freedom, the S-matrix factors do, so that the matrix element over the projectile wave function becomes an operator on the target Hilbert space. Also, it depends only on the integrated target degrees of freedomas required for our definition of the density matrix. The variables S andS in Eq. (3.9) are the analogs of j and j in Eq. (2.1).
We now want to trace the evolution by an additional rapidity ∆y so that the extra rapidity moves the observable away from the target. This can be achieved by boosting the target by an additional rapidity ∆y relative to the lab frame, or, alternatively, boosting the projectile together with the observableÔ by the same rapidity, so thatÔ remains at the fixed rapidity from the projectile. It is straightforward to do the latter. After boosting the projectile we have where for small ∆y the coherent operator C ∆y for dilute projectile (see for example Ref. [33] for the complete definition of the operator C y ) can be expanded into a power series Here b a i (x ⊥ ) is the Weizsacker-Williams field of the projectile, and j a P is the color charge density of the dilute projectile (not to be confused with j defined above which in the present context is the color charge density of the target).
The evolution equation for the operator is obtained from We remind the reader the following key identities valid for any multigluon state in the projectile Hilbert space (see Ref. [34]): (3.14) By construction, the operatorÔ commutes with the soft gluon operators a and a † , since it only involves degrees of freedom at rapidities between Y 0 + ∆y and Y + ∆y. We can thus take the averages of all the soft gluon operators in the soft gluon vacuum, since the projectile state before boost was the vacuum for these modes. Combining the expansion (3.11) and the identities (3.14) we obtain with Now recalling Eq. (3.9) we see that we can integrate the evolution kernel "by parts", so that the derivatives in the operators Q act on the target density matrix. This results in the evolution equation for matrix elements of the density matrix Note that we can rewrite this in the operator form by "integrating by parts" Q[j ] so that it does not act on ρ[j, j ] but acts on the operator whose expectation value we are calculating. Indeed let us consider the evolution of an arbitrary observable O(j, j ), Integrating by parts in the functional integral j , we arrive at Since the operator O is arbitrary, this is equivalent to the evolution of the density matrix operator in the form where the operatorQ a i [z ⊥ ] is defined in such a way that for an arbitrary ket |ψ The evolution equation for the density matrix has the celebrated Kossakowsky-Lindblad form (see the original papers in Ref. [35,36] ,ρ being the Lindbladian of the system with the so-called jump or Lindblad operatorQ a i [z ⊥ ]. The fact that the evolution has this form is not surprising. The Lindblad master equation (here without the unitary part of the evolution) is the most general form of the Markovian evolution preserving the trace and the positivity of the reduced density matrix. This equation is ubiquitous in various fields of physics whenever a description of an open system is attempted, see e.g. Ref. [37]. The meaning of the jump operator in this context is the amplitude of the process in which the "environment" experiences a quantum jump to a different level. This is very natural in the context of the high energy evolution, as Q a i [z ⊥ , j] is precisely an amplitude of emission of a soft gluon. Such emission process does indeed change the quantum state of the soft "environment".
We note one interesting feature of Eqs. (3.20, 3.17). Specifically concentrating on Eq. (3.17) we see that diagonal matrix elements ofρ evolve independently of the nondiagonal ones. This is due to the property of the operator H 3 discussed in detail in Ref. [10], valid for an arbitrary function F . Thus the diagonal matrix elements ofρ indeed evolve according to the standard JIMWLK equation.

B. Entropy growth
The Lindbladian does not have the form of the Hamiltonian evolution in quantum mechanics, since the time derivative ofρ is not given by a commutator with an Hermitian operator. Thus, the entropy ofρ increases in the course of the evolution, as well known, for the Lindblad master equation.
Here, for completeness of the discussion, we demonstrate this explicitly in the following simple way. Let us examine the effect of the evolution on a pure state and consider the evolution ofρ 2 . To reduce the notational clutter, in this section we will use the shorthand notation ρ jj = ρ[j, j ] and Q a In the first term of the equality (3.23) the k integration is trivial. Using the pure state conditionρ 2 =ρ one recognizes in this term the derivative ofρ. The rest of the terms can be rearranged after using integration by parts on the jump operators Q a k : It is clear that the evolution ofρ andρ 2 differs by a non-trivial term, indicating that a pure state becomes mixed after evolution. One can go one step further and take the trace of Eq. (3.24). Taking into account that the trace ofρ is always 1, we get Given thatρ is Hermitian and Q is real, the integrand is clearly positive definite. Therefore we conclude that d dy Thus the evolution changes the density matrixρ such that it does not correspond to a pure state anymore. In particular, by using the standard definition of the Renyi entropy S R = − ln Trρ 2 , we find d dy S R > 0, (3.27) showing that the entropy of the density matrix increases due to evolution. We will return to an explicit calculation of the Renyi entropy in the following section. The property of decoherence and (3.26) are well-known in the context of the Lindblad master equation, but has been derived in the CGC framework for the first time here.

IV. EVOLUTION IN GAUSSIAN APPROXIMATION
In the previous section we derived the evolution of the density matrix. It is highly nonlinear and non-local due to the complexity of the jump operator Q. To get some idea on how the evolution affects the off diagonal matrix elements ofρ, we will follow the ideas introduced before for the diagonal components of the density matrix in the context of the JIMWLK evolution equation, see e.g. Refs. [38,39]. Namely we will consider a Gaussian approximation for the density matrix, and will derive the evolution equations for the effective parameters. We consider the following approximation toρ: Here, to simplify the derivation, we introduced the field α instead of the color charge density. We anticipate that the rapidity evolution of the density matrix will be encoded in the rapidity dependence of parameters µ 2 y , λ y and A y . In principle µ 2 y , λ y and A y can be taken as arbitrary matrices in color space, and the form Eq. (4.1) can accommodate such a general choice. However, color neutrality requires all these matrices to be proportional to identity, and we will restrict the general ansatz correspondingly.
We start by deriving the evolution of these parameters in the dilute regime; we then also consider the approach to the saturated regime.
A. Gaussian approximation for density matrix evolution in the dilute regime To derive the evolution for the three parameters in Eq. (4.1) we have to consider three different averages Ô i and require that their evolution is reproduced by the Gaussian ansatz. The natural choice is to take the averages of the three simple linearly independent operators: .
For each operator in this set we first calculate the corresponding expectation value and then take its derivative with respect to rapidity Since the evolution of each expectation value is dictated by Eq. (3.1), Eq. (4.4) has to be equated to In fact, given our choice of operatorsÔ i it is easier to rewrite the last expression using the cyclic property of the trace as Thus the equations that determine the evolution of the parameters of the density matrix in the Gaussian approximation are Proceeding with the plan outlined above, we calculate the averages in the Gaussian state density matrix: where we have explicitly assumed that x 1⊥ = x 2⊥ . One has to be more careful with the derivation for x 1⊥ = x 2⊥ . In the dilute limit, operatorQ can be expanded to leading order in α: with corrections of order α 2 . In this limit, we also find and, finally, In what follows, in order to simplify the derivation, we set A = 0. As an initial condition, A = 0 is preserved by the evolution and is thus setting A = 0 at any rapidity is compatible with Eq. (4.14). As we will discuss in the next subsection, the function A in general does not affect the eigenvalues of the density matrix and thus can be set equal to zero for the purpose of evaluating the entanglement entropy. For a general operator, the evolution of the function A, however, can be important and should not be omitted. Note that the evolution for µ 2 y and λ 2 y decouple in the Gaussian approximation with A = 0, and the evolution equations become: These equations can be further simplified. Concentrating first on Eq. (4.15), we note that it is convenient to definē The evolution equation for this quantity becomes Physicallyμ 2 differs from µ 2 only by terms that do not depend on one of the coordinates. In momentum spaceμ 2 y and µ 2 y are therefore identical except possibly for zero momentum modes. For this reason we will not distinguish between µ 2 y and µ 2 y in the following. As for Eq. (4.16) we note that the last two terms in this equation are proportional to zero momentum modes of λ −2 y . Thus if d 2 x ⊥ λ −2 y (x ⊥ , y ⊥ ) = 0 these terms drop out. It is also easily verified that this condition is preserved by Eq. (4.16), i.e. assuming d 2 x ⊥ λ −2 y0 (x ⊥ , y ⊥ ) = 0 at initial rapidity y 0 one has d 2 x ⊥ λ −2 y (x ⊥ , y ⊥ ) = 0 at any rapidity y. We will thus drop these terms and simplify Eq. (4.16) to We observe that both Eq.  We already showed that the density matrix describing a pure initial state decoheres with evolution. Now our goal is to understand how this decoherence happens at high energies. The measure of such decoherence is the entanglement entropy. For the Gaussian density matrix, it can be calculated. Let us start with the N -th Renyi entropy, which is somewhat easier to calculate than the von Neumann one.
As we alluded to before the parameter A does not enter to the expression for the entropy. To prove this we note that in the definition of the density matrix it appears as part of a unitary basis change. In particular the density matrixρ of Eq. (4.1) can be written asρ where Thus A does not affect the eigenvalues ofρ, and does not change the evolution of any operator of the form trρ n . For that reason, in the following we will set A to zero. Using the parametrization of the density matrix we can find the N -th Renyi entropy following the same steps as in Ref. [12]. The trace of the density matrix to the N -th power is with periodic boundary conditions in the replica space α N +1 = α 1 . This integral is not diagonal in α. The easiest way to proceed with integration over replicas is to transform the expression in the exponential into the Fourier replica space. It is introduced according to This transformation satisfies the periodicity conditions α j+N = α j . The reality of α j (x ⊥ ) also leads to the relatioñ Let us consider two types of expressions that we encounter in the calculation: Using these we get The integration overα N does not involve any factors of λ. This is an integration with respect to the center of mass in the replica space. This integral will be canceled by the normalization of the density matrix. Thus we have Using the identity we obtain where T N are the Chebyshev polynomials T N (x) = cosh (N acosh x) . Thus the N -th Renyi entropy is Note that this expression gives zero in the limit λ 2 y → µ 2 y ; this can be easily established based on the property of the leading order coefficients of the Chebyshev polynomial of order N : The second term in this expansion allows to extract the first non-trivial contribution to S N close to the pure state limit λ 2 y → µ 2 y : (4.36) In the opposite limit, in a strongly mixed state, λ 2 y µ 2 y , we get The same result can be obtained from the general Eq. (4.35) taking N = 2. Here, as before, λ 2 and µ 2 are considered as operators on the transverse space.
It is clear from the above expressions that the increase of µ 2 and/or decrease of λ 2 increases the Renyi entropy and thus signals an increased mixing of the density matrix. As we have seen in the previous subsection, this is indeed what happens in the dilute regime. Using Eq. (4.20) we find that, in the dilute regime, That is, the Renyi entropy grows linearly with rapidity. Now consider the von Neumann entropy of the reduced density matrix; it is defined by S e = −Tr (ρ lnρ) .
and, assuming that lim →1 S exists, we obtain This is an exact expression for the von Neumann entropy of the reduced density matrix in the Gaussian approximation.
To demonstrate the behavior of the entropy, in Fig. 1, we plotted the function S e (x = µ 2 y λ −2 y ), where we treat µ 2 y and λ y as scalar numbers. For a strongly mixed state (i.e. |µ 2 λ −2 | 1), the evolution reads: In full generality, the evolution of the von Neumann entropy reads  Therefore, the entropy deviates fast from the pure state regime due to the presence of the logarithmic singularity in the derivative close to the pure state limit. This shows that, at least in the Gaussian approximation, a pure state quickly morphs into a mixed state.

C. Approach to saturation and Levin-Tuchin law
We now wish to study the behavior of the entanglement entropy in the saturated regime. We will again take a reasonable ansatz for the density matrix and will require that it reproduces the evolution of two simple operators. The two operators that we choose are the dipole amplitude in the fundamental representation (4.47) and the correlator As explained in Ref. [42], the operator P † (at least in the first approximation) plays the role of the conjugate Pomeron within the Pomeron field theory approximation to JIMWLK evolution. Before restricting ourselves to a particular form ofρ, let us derive the operator evolution of the two simple operators in question.

Evolution
It is straightforward to derive the evolution for operators without making simplifying assumptions about the strength of the gluon fields but instead using the full expression for the jump operatorQ a (z ⊥ ). It is in fact obvious that the dipole evolves according to the BK equation [43,44] The reason is that the operator d depends only on the eikonal matrices S, and thus its average and evolution is governed entirely by the diagonal elements ofρ in the S-basis. Since these elements evolve according to the original JIMWLK equation, so does the operator d.
The explicit calculation for P † yields To simplify this, we use S T S = 1, (T a ) T = −T a , J L S = J R and T a T a = N c , getting Interestingly we find that even in the saturated regime the charge density correlator evolves according to the BFKL equation. This is perhaps not completely surprising for the following reason. As discussed in the literature, e.g. Ref. [45], high energy evolution has a self dual structure. As a result of this dense-dilute duality, the operators that depend on the eikonal matrix S probe the structure of the target state, while those that depend on the charge operators J effectively probe the structure of the projectile state. The JIMWLK evolution describes a situation where the target is dense, but the projectile is dilute. Thus the rapidity evolution of the charge correlators reflects the rapidity evolution of charge densities in the dilute projectile, which have to evolve according to the BFKL equation.
This indeed is what we find. Note, however, that even though this result is natural, it by no means trivial, as we were only able to obtain it explicitly by using the density matrix formulation of high energy evolution, as the knowledge of the diagonal matrix elements ofρ alone is not sufficient to calculate the evolution of any function of J's.

D. The ansatz forρ in the saturated regime and the operator averages
To study the density matrix close to the saturated regime we will take a natural generalization of the Gaussian ansatz (S andS below are matrices in the fundamental representation), 52) and repeat the procedure that we performed in the dilute regime. Note that we already made the unitary transformation to rotate out the function A, as in the dilute case. Thus, similarly to the dilute regime, we have to consider only two operators in order to derive the evolution of parametersλ andμ. For simplicity, we will adopt the following natural assumptionsμ 2 . This prevents the appearance of the odderon which is not critical for our consideration relevant for high energy. In general, one can lift this assumption and repeat the derivation; for the purpose of this paper it is not necessary.
Our goal is now to calculate the averages of d and P † in this density matrix. This is not an easy task and we do not know how to perform this calculation in full generality. However, since we are interested in the behavior close to the saturation limit, we can invoke the factorized approximation used in Ref. [46][47][48]. This amounts to forgetting about the complicated group measure while integrating over S, and using the standard measure on complex numbers C for each matrix element of S. In this approximation, the averages of products of S matrices factorize into products of color singlet pairs, see Ref. [46][47][48] for details. This approximation is justified in particular when one is interested in leading powers of the area of the projectile as explained in the above references. We will not further justify this approximation here, but will instead hope that it gives qualitatively correct answers to the questions that we are asking.
Our normalization ofμ 2 is such that This means that the natural magnitude isμ 2 ∼ 1/N c . To calculate P † we will use the identities which can be trivially proven based on the definition of J a R (x ⊥ ), see Eq. (3.3). Using these identities we obtain Before acting with the second operator J a R , we note that to calculate the average we will have to setS = S after the differentiation. We will therefore only keep terms that do not vanish forS = S: Now recall that the factorization rules dictate For simplicity we calculate the averages to leading order in 1/N c . We then find which suggests thatλ −2 y is of order 1. In order to restore the natural N c power counting, we rescaleμ 2 →μ 2 Nc and obtain Recall that in the saturation regime the behavior of the dipole is governed by the Levin-Tuchin (LT) formula [49] d where ξ is a constant of order unity, and Q s is the saturation momentum. We thus conclude that in this regimē Given that the color density correlator satisfies the BFKL equation, we find whereλ 0 is determined by the initial condition. The dependence of Q 2 s on y is well known, with leading exponential behavior being where χ is the BFKL kernel, χ(γ) = 2ψ(1) − ψ(γ) − ψ(1 − γ), and γ c is the solution of equation χ(γ c ) = γ c χ (γ c ), which numerically is γ c ≈ .628. Numerically β ≈ 4.88 αsNc π [38]. Note that the density matrix is normalizable within our approximation only as long as |λ −2 | > Nc 4 |μ −2 |. Thus the calculation can only be valid for large rapidities, i.e. Since at large N c we have α s ∼ 1/N c , parametrically this is the same as Eq. (4.65) if the initial conditionλ 0 is of the order of the 't Hooft coupling,λ 0 ∼ α s N c .

E. Entropy in the saturated regime
The next natural question is how does entropy evolve in the LT regime. We can in fact adopt the results of the previous section to calculate entropy. Our approximation of calculating the functional integral corresponds simply to treating the matrix elements of S as independent degrees of freedom. The density matrix therefore behaves as a Gaussian in these degrees of freedom and the entropy is simply the entropy of a Gaussian density matrix. We therefore can directly write (4.67) Here the matricesμ 2 andλ 2 are Hermitian N 2 c by N 2 c matrices. To estimate this we need to calculate the operator Thus, interestingly the entropy grows slower in the saturated regime. This is natural since due to saturation effects the emission of soft gluons is suppressed and, thus, one expects the rate of decoherence of the density matrix to slow down.

V. WIGNER FUNCTIONAL
Let us now return to our original motivation for introducing the density matrix: can we get information on the distribution of currents in the hadronic state at high energy and, more interestingly, on correlations between currents and color charge densities? This type of question is particularly pertinent as we are interested in rare configurations in the wave function, e.g. such that produce higher than average multiplicity final states in p-p collisions. Such configurations are quite likely to also harbor large currents and therefore momentum distributions that significantly differ from the average.
A similar question in a single particle quantum mechanics is answered, at least partially, by the Wigner function, which can be approximately interpreted as giving the joint probability for the distribution of position and momentum of the particle, The momentum is proportional to the velocity of the particle, and thus the Wigner function carries information not just about the distribution of position but also about its time derivative. This joint distribution is a very interesting quantity since it can, among other things, tell us how fast the particle escapes from a given point in space.
One can define an analog of the Wigner function for a field theory. Formally let us define the Wigner functional as The high energy evolution of this functional can be readily derived from the evolution ofρ. We will not pursue this trivial derivation here but instead concentrate on a possible phenomenological application of the functional. Does this functional give us any information about the distribution of color current densities, as opposed to just the distribution of color charge densities? The color current density operator is not directly present in the effective description furnished byρ or W in Eq. (5.2). However, as it is usually the case with effective theories, certain fundamental operators can be related to objects appearing in the effective description. The only object independent of j that appears in Eq. (5.2) is the phase Φ. Thus we need to understand if Φ is related to the current density.
Let us calculate the commutation relations between the color current and the color charge density in the fundamental description. Recall that on the "microscopic" level we have the color current density [59] Commuting this with the color charge density where we used the canonical commutation relations and the Jacobi identity for the structure constants. Given this commutation relation and Eq. (5.3), we can construct operators in the effective theory which satisfy the same algebra. In particular, to reproduce Eq. (5.6) we can adopt the following representation for the current density in the "effective" description: In other words, indeed if we have a joint probability distribution of j and Φ, we also know the joint probability distribution of j and j i .
The color charge satisfies a covariant conservation equation. Disregarding for the moment the word "covariant", this amounts to This equation, in principle, tells us how fast the charge density "runs away" from any given configuration. This runaway speed is of course proportional to 1/p + , which is small. But this overall scaling is simply the consequence of Lorentz time dilation in the boosted frame. As a simple example of possible utility of the Wigner functional let us do the following simple exercise. We take a Gaussian ansatz for the density matrix Eq. (2.7), and calculate the correlation between the color charge and the color current densities. We take A = 0 for now, and obtain the Wigner functional: The normalization constant N can be obtained from the condition As one could have expected, the Wigner functional is Gaussian for the Gaussian density matrix. Using this result we can study different correlators between color currents and color densities. Below we will consider two examples. First, we start from correlators of the currents at two different positions, This correlator in a way is a proxy for two gluon azimuthal anisotropy harmonics v 2n in the CGC wave function.
Since it is proportional to λ −2 , this demonstrates the importance of the off-diagonal components of the full density matrix in relation to the momentum distribution of particles. Another illuminating example potentially pertinent to phenomenology is the correlator Eq. (5.13) demonstrates the presence of a nontrivial correlation between a proxy for the gluon multiplicity and the azimuthal anisotropy even in a simple density matrix. Again, to establish this correlation the computation of the full density matrix and not just its diagonal part was required. One feature of Eq. (5.13) is particularly interesting. Note that the correlated (connected) part of the correlator has the same energy dependence as the disconnected piece. Thus this type of correlation, if present in the wave function at initial energy, is not washed away by energy evolution.
Another interesting point is the role of the parameter A. If we reinstate it in the general Gaussian ansatz, we obtain for the Wigner functional where we introduced Although the presence of A has no effect on the correlators involving only j, it does affect j i . This is entirely analogous to how a coordinate dependent phase of a wave function does not affect the probability distribution of coordinates, but has a strong effect on the distribution of momenta (velocities) of a particle. Whether this effect is significant at high energy is an interesting question worth exploring. The upshot of this short discussion is that the Wigner functional does have a potential of being a useful tool in understanding the dynamical structure of the hadronic wave function. The quantitative study of the properties and the evolution of the Wigner functional deserves a serious effort, and is left for future work.

VI. DISCUSSION
In this paper we have introduced the notion of the CGC density matrixρ. This is the reduced density matrix in the CGC effective theory obtained by tracing over all the degrees of freedom in the QCD Hilbert space except the rapidity integrated color charge density. We stress again that this is not the same density matrix as considered in Refs. [12,13], where the valence degrees of freedom were integrated out to obtain the reduced density matrix on the soft gluon Hilbert space.
We have derived the evolution equation for the density matrixρ and have shown that it is of the Kossakowsky-Lindblad form with the jump operator being equal to the single soft gluon production amplitude. This is intuitively quite agreeable, since in general the meaning of the jump operator is to introduce a jump to a different quantum state of the "environment", which in our case contains soft gluon degrees of freedom.
The Kossakowski-Lindblad form is the most general form of Markovian evolution allowed by a probabilistic interpretation of the density matrix, i.e. overall normalization and positivity of all eigenvalues. This suggests that the general form of the evolution equation should persist beyond leading order. It is thus possible that one can simplify the derivation of the NLO JIMWLK [50][51][52][53] by directly calculating corrections to the jump operator, rather than to the JIMWLK Hamiltonian, which is a more complicated object. One may hope that the same framework can also accommodate improved leading order JIMWLK versions which resum large transverse logarithms. Physically one expects that since the evolution in energy is aligned with the evolution in the frequency of produced gluons, the typical time scale of the evolution will remain always larger than the time scale of the soft gluon fluctuations. If this is the case, the Markovian nature of the evolution should be preserved beyond the leading order. Although physically reasonable, a better understanding of possible sources for non-Markovian effects in the evolution is necessary.
The Kossakowski-Lindblad evolution is known to lead to increasing entanglement entropy with the evolution "time". We have indeed calculated the evolution of entanglement entropy in a Gaussian approximation, both in the dilute regime and close to saturation. We found that in both cases the entanglement entropy increases linearly with rapidity. In the dilute regime the rate of increase coincides with the leading BFKL eigenvalue, while in the dense (Levin-Tuchin) regime it is half of that value. The slower growth of entropy in the saturated regime is likely caused by the suppressed emission probability of soft gluons close to saturation.
The linear growth of entropy with rapidity is a rather interesting result. One may naively expect that the entropy associated withρ is proportional to the total gluon number n y -at least as long as n y is not too large. However this is not the case. The total number of gluons in the dilute regime grows with rapidity exponentially, while the entropy Eq. (4.44) only grows linearly and thus much more slowly. The same type of behavior persists in the dense regime.
A similar behavior of the entanglement entropy of the proton in the context of DIS was proposed in Ref. [54]. The picture of Ref. [54] is very simple. It assumes that all partonic states in the proton wave function at high energy completely decohere from each other and all accessible states become equally probable. Thus the density matrix becomes proportional to the unit matrix on the subspace of the Hilbert space which is "populated" at a given energy. The dimension of this subspace is proportional to the mean number of gluons in the wave function, which grows with the BFKL exponential, d ∝ e γy . The normalization ofρ means that it has e γy equal eigenvalues, each one approximately ρ i ∝ d −1 . For such a density matrix we know that the entropy S e ≈ − ln ρ i = γy. The growth of this entropy with energy is slow becauseρ is already maximally mixed on the subspace of dimension d, and the growth is only due to the increase of this dimension with energy.
Although we do not know how closely the density matrix introduced in the present paper is related to the object considered in Ref. [54], it is instructive to examine our formulae and understand whether the behavior we find indeed conforms with this simple argument.
Consider for simplicity the density matrix Eq. (4.23). Although above we have used this Gaussian ansatz only in the dilute regime, the following qualitative discussion should apply to both regimes. At very high energy, i.e. close to saturation, λ −2 y is very large and the density matrix Eq. (4.1) indeed becomes very close diagonal ρ(α, α ) ∝ δ(α − α ). Let us for the moment assume that we indeed can neglect the nondiagonal matrix elements ofρ.
Then since µ 2 y is also very large, for values of the field α such that α 2 < µ 2 y the matrix elements ofρ do not depend on α. Thus, in the high energy regimeρ in effect is proportional to a unit matrix of dimension d ∝ |α max | ∝ e γ 2 y . The entropy associated with such a density matrix should be given by ln d = γ 2 y. Although qualitatively correct, we are missing here a factor of 1/2 relative to our result Eq. (4.44). A closer look at our derivation indeed reveals the origin of the missing factor 1/2. As is obvious from Eq. (4.43), only half of the entropy growth comes from the growth of µ 2 and therefore of the dimension d. The other half is contributed by the increase of λ −2 , which controls the extent to which off diagonal elements ofρ are negligible. Therefore, in the dilute regime (but at high enough energy where e γy 1) the entropy grows due to two distinct effects: growth of the dimension d of the subspace on whichρ is nonvanishing, as well as further decoherence ofρ on this subspace. The two effects contribute equally to the entropy.
We note that we did obtain S e ≈ γ 2 y in the LT regime, and one might think that in this saturated regime the previous argument holds. However, a closer inspection shows that this is not the case. In the saturation regime, just like before the dimension of the relevant Hilbert space on which the density matrix is close to unity is controlled by the parameter µ 2 . However µ 2 now grows with energy much slower than exponentially, i.e. Eq. (4.61). Thus, the "expansion" of the populated subspace of the Hilbert space does not contribute any linear in rapidity term to the entropy. On the other hand, the growth of λ −2 is still exponential just like in the BFKL regime, Eq. (4.63). All the entropy evolution in Eq. (4.70) therefore originates from further decoherence on the Hilbert space of approximately fixed dimension. Therefore, neither in the dilute nor in the dense regime, the linear growth of entropy discussed in the present paper seems to originate entirely from a picture proposed in Ref. [54], although this conclusion could be basis dependent.
As the last point in this paper, we have also defined the Wigner functional associated with the density matrixρ. We have argued that it can give access to understanding the joint probability distribution in the space of color charge density and color current density. We hope that such a distribution can teach us about the momentum distribution of produced gluons in events with high multiplicity, which should be instrumental in understanding the correlated behavior of produced hadrons. We have shown that even within the simple Gaussian ansatz, the distributions of color charges and color currents do not factorize, and that this non factorization feature is not eliminated by high energy evolution.
We hope that further work on these subjects will lead to a better, and more complete understanding of hadronic physics at high energies.