JIMWLK evolution, Lindblad equation and quantum-classical correspondence

In the Color Glass Condensate (CGC) effective theory, the physics of valence gluons with large longitudinal momentum is reflected in the distribution of color charges in the transverse plane. Averaging over the valence degrees of freedom is effected by integrating over classical color charges with some quasi probability weight functional W [j] whose evolution with rapidity is governed by the JIMWLK equation. In this paper, we reformulate this setup in terms of effective quantum field theory on valence Hilbert space governed by the reduced density matrix ρ̂\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \hat{\rho} $$\end{document} for hard gluons, which is obtained after properly integrating out the soft gluon “environment”. We show that the evolution of this density matrix with rapidity in the dense and dilute limits has the form of Lindblad equation. The quasi probability distribution (weight) functional W is directly related to the reduced density matrix ρ̂\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \hat{\rho} $$\end{document} through the generalization of the Wigner-Weyl quantum-classical correspondence, which reformulates quantum dynamics on Hilbert space in terms of classical dynamics on the phase space. In the present case the phase space is non Abelian and is spanned by the components of transverse color charge density j. The same correspondence maps the Lindblad equation for ρ̂\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \hat{\rho} $$\end{document} into the JIMWLK evolution equation for W .


JHEP05(2020)036
3 The reduced CGC density matrix and its evolution 8 3.1 Defining the charge shift operator 9 3.2 The evolution 11 4 The differential form of the evolution -the Lindblad equation  34 5.4 From JIMWLK back to Lindblad 37 1 Introduction This paper examines the status of the JIMWLK evolution equation [1][2][3][4][5][6][7][8] in relation to the effective density matrix of a high energy hadronic system. We are motivated to address this question by the discussion in a recent paper [9] which suggested an extension of JIMWLK framework to include a wider set of observables other than just color charge density j a (x) in the hadronic wave function. The starting point of [9] is the interpretation of JIMWLK evolution equation as the equation for diagonal matrix elements of the density matrix in the JHEP05(2020)036 color charge density basis. Although this interpretation is natural when the color charge density is large, it is not quite clear how to formalize it for low density, since in this regime the commutator of the color charge density operators is not negligible and a basis in which all components of j a (x) are diagonal obviously does not exist. On the other hand, as shown a while ago [10] the calculations of averages in this regime as well can be formulated in terms of the functional integral over classical fields j a (x), which suggests that perhaps such interpretation albeit possibly modified, can be put forward after all. An interesting suggestion of [9] is that the rapidity evolution of the generalized CGC density matrix is of the Lindblad type [11,12]. This type of evolution is very general in quantum mechanical systems where one follows only part of the degrees of freedom by reducing the density matrix over the "environment" (the unobserved degrees of freedom in the Hilbert space). If the "environment" degrees of freedom have only short range correlation in time, the dynamics of the observed part of the system is Markovian and is therefore governed by a differential equation. The Lindblad form of such evolution is the most general one that preserves the properties of the density matrix stemming from its probabilistic nature (normalization and positivity of all eigenvalues). Although Lindblad equation naturally arises in time evolution of quantum systems, JIMWLK evolution is of a somewhat different nature. It describes the change of the system with rapidity (or energy) but not in time. It is thus not obvious whether one should expect Lindblad form to be generic in this context and if yes, under what conditions.
In this paper we try to address these questions. We arrive at two basic results. First, we show that JIMWLK evolution can indeed be understood as evolution of a density matrix. Within the JIMWLK framework however, the density matrix is not generic, but is rather assumed to depend only on the operatorsĵ a (x) which satisfy the standard SU(N ) commutation relations. The fact thatρ depends only on the generators of the SU(N ) group means that it has a quasi diagonal form -i.e. it does not connect states belonging to different representations of SU(N ). It is in this sense that the reduced density matrix is (almost) diagonal even if the commutators of j a cannot be neglected. The consequence of this strong assumption on the nature of the density matrix is that the states in the Hilbert space of the reduced system are completely specified by their SU(N ) transformation properties at every transverse position x, and therefore in the reduced space one loses track of the gluon longitudinal momentum as well as polarization.
Second, we show that on this Hilbert space the JIMWLK evolution is indeed equivalent to Lindblad type equation for this restricted set of density matrices. The same applies to the so called KLWMIJ evolution which describes the dilute regime. In analogy with time evolution of quantum mechanical systems, the Lindblad equation arises in the situation when the correlations in the "unobserved" part of the system are short range in rapidity. However we also argue that in general (i.e. away from the dense -JIMWLK and dilute -KLWMIJ limits) the high energy evolution is unlikely to be of Lindblad type. This follows from certain general properties of our derivation of the evolution of the density matrix based on the calculation of the CGC wave function presented in [13]. Although the calculations of [13] are strictly valid only in the aforementioned limits, the general features of the derivation are expected to be more universal. The reason that the Lindblad JHEP05(2020)036 form is not expected to arise, is that in the high energy evolution framework, the rapidity does not just play the role of the evolution parameter, but also that of the label on the quantum states of the gluons which are integrated out. In this situation in general one does not expect the Lindblad form for the differential equation. Thus to ensure Lindblad form nontrivial conditions on dependence of the matrix elements on gluon rapidities have to be satisfied. We discuss this point in detail in section 4.
Another result of this paper is the precise mathematical relation between the effective density matrix, and the "probability density function" W [j] that usually appears in the literature as the subject of JIMWLK evolution. We confirm that the quantum mechanical averaging with the density matrixρ is mapped into the calculation of observables in terms of functional integral over classical fields j a (x) with the weight functional W [j], as indeed always done in the CGC literature. This functional integral must be regarded as an integral over the phase space variables of the classical system, and not its configuration space variables. This quantum-classical correspondence between the quantum density matrix and the classical functional of phase space variables W [j] is deeply related to the correspondence between the density matrix and Wigner function in ordinary quantum mechanics. In the context of high energy evolution we require a generalization of the original Wigner-Weyl correspondence [14] since the phase space of the theory is spanned not by operators q and p which constitute the Heisenberg algebra, but rather by operators j a with the SU(N ) algebra. Nevertheless the basic correspondence involves mappings between quantum operators and Hilbert space on one side and classical quantities and phase space on the other side in the sense of Weyl's correspondence rule. The weight functional W [j] is consequently identified as the Wigner functional [14] and can indeed be considered as a quasi-probability distribution on the phase space.
The outline of this paper is the following. In section 2, we give a brief review of the Lindblad equation for density matrix of an open quantum system and a recap of the Hamiltonian formalism of CGC effective theory. In section 3 we explain how to define the reduced CGC density matrix, and show that its rapidity evolution has the Kraus form, which is a general evolution that preserves probabilistic interpretation of a density matrix, but is not necessarily differential. In section 4 we derive the differential evolution of the density matrix using the analog of Markovian porperty, i.e. the fact that the correlation of the "environment" is short range in rapidity. We show that in the dilute (KLWMIJ) and dense (JIMWLK) limits the differential evolution equation is indeed of the Lindblad type. We also discuss the properties of the derivation which suggest that the standard Lindblad form is bound to be modified away from these limits. To be clear, in this paper we do not go beyond the original JIMWLK setup in the sense that we consider density matrices that only depend on color charge density operators, and thus presently our derivation does not extend to ideas put forward in [9]. In section 5, we derive the explicit relation between the standard JIMWLK approach and the density matrix approach described in this paper. We show that the two are related by a variant of the Wigner-Weyl quantumclassical correspondence and spell out explicitly the correspondence rules which transform one setup into the other. The JIMWLK and KLWMIJ equations are then reproduced by mapping the Lindblad equation for the density matrix in the appropriate (dense and dilute) limits into the classical phase space. Finally section VI contains a short discussion.

Lindblad equation for open quantum systems
In this section we present a short review of Lindblad equation for open systems.
Lindblad equation is the most general Markovian and non-unitary evolution equation for density matrix of an open quantum system. This equation preserves the properties of density matrix: hermiticity, unit trace and positivity. Here we follow the heuristic discussions by Preskill [15]. More physical derivations and applications can be found in the books [16,17]. Consider a bipartite system involving two subsystems: the "observed system" and the "environment" with the HamiltoniansĤ s andĤ e , respectively. The two subsystems interact via the HamiltonianĤ se . The total density matrix of the complete system evolves according to the quantum Liouville equation withĤ =Ĥ s +Ĥ e +Ĥ se . Formally, the solution can be expressed aŝ To obtain the density matrix of the observed subsystem after a finite time, one traces over the Hilbert space of the environment. Let us assume that the initial total density matrix is a direct product of the density matrices of the observed system and the environmentρ(0) =ρ s (0) ⊗ρ e (0) =ρ s (0) ⊗ |0 e 0 e |. For simplicity let us take the environment to be initially in a pure state denoted by |0 e , which can be thought of as the ground state without loss of generality. The density matrix of the observed system is then expressed aŝ Here {|n } represents a complete basis in the Hilbert space of the environment. The objectŝ M n (t) = n|Û (t)|0 e , sometimes called superoperators, are operators on the Hilbert space of the observed system and govern the evolution of its density matrix. As far as the dynamics of the environment is considered,M n (t) represents the transition amplitude for the environment, which is initially in the state |0 e , to be in the state |n after a finite time t. They satisfy the property nM

JHEP05(2020)036
The Markovian approximation holds if the typical correlation time between the environment degrees of freedom t corr is shorter than the typical inverse frequency of the observed system ∆t s , which is of the order of the relevant "discretization" time step for approximate differential time evolution. If this is the case the state of environment is only affected by the state of the observed system at the particular time of observation (measured with accuracy ∆t s ), and thus the back reaction -the effect of the environment on the observed system is local in time. We note that this is the typical Born-Oppenheimer situation, when the environment is associated with fast degrees of freedom, while the observed system is relatively slow. In the opposite regime it is clear that local (differential) time evolution is impossible, since the backreaction of the environment on the system will depend on the state of the system at some past time.
In Markovian regime one then proceeds as follows. For an infinitesimal period of time, only terms linear in dt should be kept on the right hand side of eq. (2.3). The superoperators for n > 0, have the structureM The argument here is thatM † n (t)M n (t) is the probability for the environment to "jump" to the state n during the time t. For small enough t (but such that t > t corr ) this probability should grow linearly with t. The operatorsL n are called Lindblad operators or jump operators as they involve transitions of the environment to different states after an infinitesimal time.
The remaining superoperator has the form M 0 (dt) = 1 + (−iĤ s +K)dt (2.5) with H s and K being Hermitian. This is the transition amplitude for the environment to be in its original state after an infinitesimal time and should be linear in time for small enough times. The operatorK is related to the wave function renormalization effect andĤ s governs the unitary evolution of the system without causing any changes to the environment. The Kraus normalization condition nM † n (dt)M n (dt) = 1 relates the wavefunction renormalization operatorK to the jump operators bŷ Taking the limit dt → 0, the Kraus representation then becomes an differential equation This is the Lindblad equation, or sometimes known as Gorini-Kossakowski-Lindblad-Sudarshan master equation [11,12]. JHEP05(2020)036

The soft gluon vacuum and the CGC
We now review the derivation of the high energy evolution [18]. There exist two equivalent approaches to the derivation of the CGC effective theory. One is based on the Lagrangian formalism [2-4, 6, 7, 19, 20] and the other on the Hamiltonian formalism [21,22]. We briefly review the Hamiltonian formalism as it will be the starting point for deriving the Lindblad equation for the CGC density matrix. The derivation of the JIMWLK evolution equation starts with establishing the ground state wave function of soft gluon modes in the background of more energetic gluons which are described by a color charge density field.
In the light cone gauge A + = 0, the Hamiltonian of the pure gluonic sector of QCD is with the chromoelectric and chromomagnetic parts being The covariant derivative is defined as D ab i = ∂ i δ ab − gf acb A c i and ∂ − = ∂/∂x − is the longitudinal spatial derivative. The 1/∂ − operator in the expression of the chromoelectric field has to be regularized as it contains a singularity at vanishing longitudinal momentum, k + = 0. This singularity is ultimately related to the zero mode in the A a i (x − , x ⊥ ) fields and is regulated by imposing a residual gauge fixing condition. We choose the residual gauge fixing One separates the gluonic degrees of freedom imposing a longitudinal momentum separation scale Λ + . In the high energy limit, the dominant interaction between soft gluons (k + < Λ + ) and hard gluons (k + > Λ + ) has the form of eikonal coupling A − a J + a with J + a representing the color charge density of the hard gluons and A − a representing soft gluons. This interaction term emerges from the chromoelectric part of the Hamiltonian and involves the specific expressions Furthermore, as far as soft gluons are concerned, the hard gluon dynamics can be viewed as frozen in time so that the color current J + a ≡ J + a (x − , x ⊥ ) is time independent at the lowest order. All in all, the Hamiltonian for the soft gluonic modes becomes Canonical quantization is implemented by promoting the normal modes of the full A a i fields to operators and imposing the equal (light cone) time commutation relation

JHEP05(2020)036
with the sign function defined as (x) = 1 2 (Θ(x) − Θ(−x)). In terms of the canonical creation and annihilation operators, the normal modesÂ a i have the expansion 14) The color charges in the leading order are taken to have the extreme Lorentz contracted The components of color charge satisfy the commutation relations of the SU(N ) algebra (2.16) The Hamiltonian system eqs. (2.11), (2.10), together with the commutation relations eqs. (2.12), (2.16) constitute the starting point for the derivation of the CGC effective theory. The first goal is to find the ground state wave function of the soft glue. This is in general a very complicated problem, but is simplifies in two parametric regimes. One interesting regime is when the color charge density is small j a ∼ g (dilute limit). Here one can treat the interaction with color charges perturbatively. The other regime is the dense limit where the color charge is parametrically large j a ∼ 1 g . Here the simplification is that the commutator of the color charges is (almost) negligible and they can be treated as (almost) classical fields.
In the dilute limit, the ground state wave function can be found by a direct perturbative calculation. The resulting vacuum wave function can be written as with the coherent operator where E is the energy of the process.
Here b a i (x ⊥ ) satisfy the equations 19) or, in the dilute regime (2.20)

JHEP05(2020)036
In the dense limit, similar analysis applies except now b a i (x ⊥ ) ∼ 1/g and additional order O(1) quantum fluctuations on top of the b a i fields need to be considered. One can still use perturbative expansion in g, but resumming terms of order gb. In the leading order the Hamiltonian is diagonalized by a nontrivial Bogoliubov transformation. The detailed analysis appears in [22]. The resulting ground state wavefunction is The additional Bogoliubov operatorB can be formally expressed aŝ Here α, β represent all the possible indices (color, spatial coordinates, polarization, and longitudinal momentum which varies between Λ + e −∆y and Λ + ). The explicit expression of the symmetric matrix M αβ is not available, however, the action of the Bogoliubov operator on the fundamental degrees of freedomÂ a i andĵ a have been derived. The nontrivial structure of the soft gluon ground state leads to appearance of induced color charge density due to the soft gluons modes. This additional color charge density serves as an additional source for even softer gluons which arise in the evolution to even higher rapidities. This is the basic physics of the high energy evolution.

The reduced CGC density matrix and its evolution
Having found the vacuum of the soft gluons, we can now address the evolution at high energy. We take here a different perspective on this derivation than given in the literature, and discuss the evolution from the point of view of quantum density matrix.
Given that we have separated our degrees of freedom into soft and hard gluons, we can view our system naturally as bipartite. At some initial rapidity, the soft gluons are in the perturbative vacuum state, and thus the total density matrix is separablê where the density matrixρ v is an operator on the hard gluon Hilbert space.
The assumption inherent in the derivation of the JIMWLK equation is that the only relevant degrees of freedom on this Hilbert space are components of the color charge densitŷ j a (x ⊥ ). This is a crucial assumption. If the valence Hilbert space could be factorized into a direct product of the space spanned byĵ a (x ⊥ ) and its complement, reducing over the complement would rigorously defineρ v [j]. However the full Hilbert space of the valence modes does not have such a direct product structure. It is thus not clear whether a well defined mathematical procedure of "integrating out" exists which may reduce the density matrix so that in general it depends only onĵ a (x ⊥ ). Nevertheless one can simply assume that at initial rapidity the density matrix indeed has such a form. It is then true (as we will see below) that this form persists throughout the evolution to higher rapidities. We will thus abide by this assumption and will treatρ v as an operator that depends only onĵ a (x ⊥ ).

JHEP05(2020)036
After boosting the system by a finite rapidity ∆y, the total density matrix changes due to the emission of soft gluons into the newly opened rapidity interval.
The gluon emission operator as discussed above can be written aŝ withĈ andB defined in eqs. (2.18) and (2.22), respectively. This form applies both for the dilute and dense regime of the evolution. Note thatΩ depends on the soft gluon creation and annihilation operators as well as the color charge density operator. While theâ a † i ,â a i act on the soft vacuum state |0 , theĵ a acts on the valence (hard) density matrixρ v . Dependence on ∆y ofΩ is crucial in obtaining the evolution equation. This point will be elaborated in the following. Our next goal is to derive the reduced density matrix by tracing over the "environment" degrees of freedom. The purpose of this reduction of the Hilbert space is to integrate out all the additional degrees of freedom associated with soft gluons that emerged after boosting the wave function, except the additional color charge density that they contribute. The reason for this exception is, that in the next step in the evolution the even softer gluons will couple to the total color charge density, including that due to gluons in the rapidity interval between y and y + ∆y. Our current soft gluons give a nontrivial contribution to this charge density, and we have to keep this extra contribution explicitly, rather than integrate it out. Hereĵ a soft (x ⊥ ) has the explicit expression eq. (2.15) with the longitudinal momentum integration restricted in the rapidity range ∆y.
It is thus clear that we should not simply reduce the density matrix over the Hilbert space of soft gluons, but "partially" trace over the soft gluons integrating out all degrees of freedom except the color charge density. To facilitate this partial tracing over soft gluons, we introduce the operatorR, which is defined by its action onĵ a , We will look forΦ a (we omit the transverse coordinate dependence for simplicity) as a set of operators acting on the same Hilbert space asĵ a satisfying the following commutation relations with M chosen to satisfy the requirement exp iĵ a softΦ a ĵ e exp −iĵ a softΦ a =ĵ e +ĵ e soft . (3.10) In calculating the action ofR we assume that the operatorsĵ a satisfy SU(N ) algebra, and so do the operatorsĵ a soft , while the two set of operators commute with each other. We use the Baker-Hausdorff formula With the commutation relations eq. (3.9) we have (for adjoint representation −if abc = T a bc ) (3.14) Let us take the ansatz Clearly, c 0 = 1 follows from the requirement eq. (3.10). This requirement further imposes the constraint which after substituting the ansatz for M (χ) becomes

JHEP05(2020)036
which is equivalent to the following recursive relations These relations are satisfied by One can check explicitly that Taylor expansion of eq. (3.19) in χ reproduces all the coefficients calculated using the recursive relations in eq. (3.18).
Once the function M in eq. (3.19) has been determined, the algebra ofĵ a andΦ a is completely defined.
We note that in order for this algebra to be consistent, the commutators have to satisfy Jacoby identity 20) which is equivalent to an additional constraint on M In appendix A we verify that the Jacoby identity is in fact satisfied at least to fourth order in expansion of eq. (3.19) in powers ofΦ. Although we do not have a complete all order proof, one could in principle continue the order-by-order proof. We believe that the algebra eqs. (3.9), (3.19) is in fact consistent and we will continue our analysis under this assumption. We have thus found the algebra of operatorsΦ a andĵ a that implements eq. (3.7). Note that expansion of M in powers ofΦ can be recast as formal expansion ofΦ in powers of δ δĵ a . Thus to leading order we haveΦ a = −i δ δĵ a + . . .. In the dense regime where the commutators of charge densities can be neglected, our operatorR therefore reduces precisely to the shift operator exp{− x ⊥ j a soft (x ⊥ ) δ δj a (x ⊥ ) } used extensively in the existing literature. The previous discussion puts its use also in the dilute regime on firm mathematical basis, provided the commutation relations ofΦ a are modified according to eq. (3.19).

The evolution
Having defined the charge density shift operatorR we can write eq. (3.5) in the form The operatorR in this expression can be understood as acting on the density matrixρ v rather than on the observable O. Using this form we can define the reduced density matrix which when traced with the operator O(ĵ) gives the same result asρ(∆y) traced with

JHEP05(2020)036
O(ĵ +ĵ soft ). We thus define the evolved CGC reduced density matrix by tracing over soft gluonŝ withM n = n|RΩ|0 . The complete basis {|n } represents the Fock states in the soft gluon Hilbert space. This procedure technically is very similar to the standard reduction of the Hilbert space discussed in the previous section withρ v (∆y) playing the role of the reduced density matrix in a bipartite system. When formulated in this way, the rapidity evolution of the CGC density matrix is formally very similar to the time evolution of the reduced density matrix of a bipartite system with the operatorRΩ playing the role of the time evolution operatorÛ .
Eq. (3.23) gives the change of density matrix in the form of a Kraus representation. As a consequence,ρ v (∆y) has all the properties of a density matrix as long asρ v is a density matrix initially. Note that, nM † nMn = 1 as bothΩ andR are unitary operators.

The differential form of the evolution -the Lindblad equation
To extract a differential equation from the Kraus representation, we need to evaluate the superoperatorsM n and analyze their ∆y dependence. The calculations ofM n can be simplified by working in the Leading Logarithmic Approximation (LLA) so that only terms that are proportional to α s ∆y on the right hand side of eq. (3.23) are kept.

The dilute limit
We start by considering the dilute limit, i.e. assume that parametrically b a i ∼ O(g). In this regime the gluon emission operator is just the coherent operator and Λ = 0. This is the so called KLWMIJ limit introduced in [23].
Note that we have changed the integration variable from longitudinal momentum k + to rapidity η and an explicit numerical factor √ 2 follows [13]. In this limit the dependence on ∆y becomes very transparent Note that the operatorR has a nontrivial action on the n-gluon state. It does not change the number of soft gluons in a Fock state but rotates their color indices according to its definition in eq.
In the LLA we need to collect terms which contribute at order O(α s ) to the evolution. For the virtual term we haveM It is obvious that, for Fock states with even numbers of gluons,M 2m is at least of order O(g 2 ) and thus will not contribute to the evolution at LLA. The same holds forM 2m+1 associated with Fock states of odd numbers of gluons. The only exception isM 1 related to the single gluon Fock state. For a one-gluon Fock state |1 {α 1 ,w 1 ,η 1 } = a † α 1 (w 1 , η 1 )|0 with transverse position w 1 , rapidity η 1 and color index α 1 we have, Summing over all possible one-gluon Fock states, withb α = R † αβ b β and again we used the compact notation with index α, β representing colors, transverse coordinates and polarizations. The evolution equation for the density matrix follows In this equation we have written the virtual terms in terms ofb rather than b, since the unitary operator R drops out of this expression anyway. This is the Lindblad equation for the CGC density matrix in the dilute limit. Eq. (4.9) is written in a somewhat convoluted form in terms of the operatorsb, which contain the operator R. It is perhaps worthwhile to make explicit the operational meaning of various factors of R in the right hand side of eq. (4.9). As already mentioned, the virtual terms do not actually involve R since for unitary R (4.10) As for the real term, we have (suppressing the transverse coordinate) where the last term is defined by Taylor expanding ofρ v , shifting the argumentĵ a in every term by the matrix T a and finally taking the γβ matrix element of the the whole expression in all products of T 's that arise.

JHEP05(2020)036
In the above explicit calculation, the LLA automatically picks up terms that are linear in ∆y thus making the extraction of a differential equation from the Kraus representation straightforward. Physically indeed we can understand this from the point of view of Markovian nature of the process. The variable analogous to time in the present discussion is rapidity. Thus the requirement of short range correlations in time of the "environment" in the CGC case translates into the requirement of short range in rapidity correlations for the soft gluons, which are integrated over. This is indeed the case. In the LLA the relevant "time" scale for the change of the density matrix is O(1/α s ), as obvious from the differential equation eq. (4.9). The soft gluons in our approximation do not interact with each other, and thus their correlation function is free. The free propagator is proportional to 1/k + ∼ e −y , and thus the typical correlation length in rapidity space is O(1). The evolution is therefore clearly in the Markovian regime which allows, at least naively speaking for the existence of differential evolution in the Lindblad form. We will come back to the discussion of Lindblad form later.
Equation eq. (4.9) may look slightly unfamiliar as it does not quite have the form of the KLWMIJ equation discussed in [23]. This is because it is written for density matrix and not the weight functional W [j]. To get to the latter form one needs to perform an extra step, i.e. Weyl transformation. This will be the subject of the next section. But before we do that, we consider the evolution of the density matrix in the dense regime.

The dense limit
As we have seen, in the limit where the hadronic wave function contains a small number of partons (the dilute limit), the Lindblad form of the evolution equation follows directly using the straightforward perturbation theory at low x. We now turn our attention to the dense limit, where we assume that the color charge density in the wave function is large, parametrically of order 1/g. The wave function in this limit has been calculated several years ago in [13]. In this section we use the results of that paper and reinterpret them from our current point of view.
To prepare for the calculation, note that the soft gluon emission operatorΩ, when acting on the vacuum state |0 can be written aŝ Here we write out the dependence on rapidity explicitly. Other indices (color, polarization, transverse position) are collectively represented by the Greek letters α, β, γ. The matrix Λ ab ij (x ⊥ , y ⊥ , η 1 , η 2 ) determines the amount of "squeezing" of the soft gluon vacuum. As we mentioned above, it has not been calculated explicitly in [13], however its properties relevant to the JIMWLK limit are known (see later). The N (Λ) is a normalization constant that only depends on Λ. Note that both b α and Λ βγ are operators in the Hilbert space of hard gluons as they depend on the color charge density j, and so in principle they do not commute. In eq. (4.12) all the factors of Λ should be understood as placed to the right JHEP05(2020)036 of b α . In the JIMWLK limit however, where parametrically, b = O(1/g) while Λ = O(1), as was shown in [13] the order of the factors does not matter. In fact in showing that the operatorΩ eq. (4.12) diagonalizes the QCD Hamiltonian to leading order, refs. [13,22] explicitly used this argument and assumed commutativity of the various factors of b and Λ. We will not deviate from this assumption here and will treat these factors as commuting.
We further separate the annihilation operatorâ α (η) from the coherent state operator and move it to the far right acting on the vacuum state: where we have defined (4.14) Since Λ αβ (ζ, η) depends only on the rapidity difference ζ − η [13], Λ 0,αβ is rapidity independent. It does however have a nontrivial dependence on the width of the evolution step ∆y. The nature of this dependence is very important. As we discussed above, we expect to have a bona fide differential evolution equation only if the correlations of the soft gluons in rapidity are short range. The function Λ(η, ξ) is in fact the inverse of the correlator of the soft gluon modes. It should therefore decrease exponentially for rapidity difference greater than ∼ 1. For such a function Λ the dependence of Λ 0 on ∆y should be smooth with Λ 0 approximately constant for 1 < ∆y < 1/α s . We will assume here that this is indeed the case and will treat Λ 0 as a constant independent of ∆y. The results of [13] suggest that this is valid in the JIMWLK limit, i.e. when the dense hadron scatters on a dilute target, which is the regime that concerns us in this paper. We note that going beyond the JIMWLK limit posed some problems in [13], precisely for the reason that some of the soft modes in general seemed to possess long range correlations in rapidity. Our current understanding is that such long range correlations indeed are incompatible with the differential form of the evolution. It is thus possible that in order to go beyond the JIMWLK limit one would have to rethink the way in which the bipartitioning into the "observable" system and "environment" is done. This is however beyond the scope of the present paper. In eq. (4.13), the first exponential represents wavefunction renormalization effects that have an overall ∆y factor. The second exponential contains the single gluon emission vertex ib α (1 − Λ 0 ) αβ which is "renormalized" relative to the dilute case by the presence of the Bogoliubov operator B, while the third exponential contains the double gluon emission vertex Λ αβ (ξ, η).
Two gluons emitted from the same double gluon emission vertex are in general correlated in rapidity, while two gluons emitted from two single gluon emission vertexes are uncorrelated.
We are now ready to calculate the superoperators. The fundamental difference with the dilute case, is that now not only one gluon state, but states with arbitrary number of JHEP05(2020)036 soft gluons yield nontrivial jump operators that contribute to the evolution of the density matrix. For an n soft gluon state we havê Depending on the Fock state |n being considered, we separately discuss the situations when the Fock state contains zero gluons, odd number of gluons and even number of gluons.

Wavefunction renormalization operator
The superoperatorM 0 represents the wavefunction renormalization effectŝ Up to terms linear in ∆y, Note that the wavefunction renormalization operatorM 0 is independent ofR and we have ignored the normalization N (Λ) factor, since it is irrelevant in the JIMWLK limit [13]. The superoperatorM 0 contributes to the change of density matrix through the term

Jump operators with odd number of gluons
For Fock states with odd numbers of gluons, one needs odd number of single-gluon-emission vertices in calculating the jump operators. However, every single gluon emission brings an extra power of ∆y, since gluons produced from different single-gluon-emission vertices are uncorrelated in rapidity. The integral over rapidity of every such gluons in the amplitude and conjugate amplitude brings therefore an extra power of ∆y. Thus one needs to keep only one single-gluon-emission vertex inM 2i+1 in order to calculate the relevant jump operators that contribute to the differential form of the evolution equation. The explicit expression for a jump operator follows from eq. (4.15) HereΛ

JHEP05(2020)036
To arrive at this expression we have inserted the factorR †R = 1 next to the soft gluon vacuum state |0 , used the fact thatR|0 = |0 and evaluated the action ofR on the soft gluon creation and annihilation operators using eq. (4.3).
Importantly, the operator ordering in eq. (4. 19) is such that all the operators R are understood to be placed to the left of all the factors of the b α and Λ αβ . This follows from the fact that the operatorR in the original expression is acting directly on the n-gluon state, and thus all the factors ofΦ indeed are ordered to the left of all j-dependent factors in the original expression. Thus for example in the definition eq. (4.20) the action of R δβ on Λ is understood only as a color matrix rotation. This comment also applies to the rest of the formulae in this section.
In the following, we explicitly calculate a few expressions of the jump operators and their action on the density matrix. This will make the dependence on ∆y more transparent.
For a one-gluon Fock state |1 {α 1 ,w 1 ,η 1 } = a † α 1 (w 1 , η 1 )|0 with transverse position w 1 , rapidity η 1 and color index α 1 , the jump operator iŝ Note that the jump operator associated with one-gluon Fock state is independent of the rapidity index η 1 . Integration over all the one-gluon Fock states produces an overall factor ∆y in the evolution of the density matrix. The one-gluon jump operators contribute to this evolution througĥ In the last line we have reverted to the convoluted notation where single index α represents the transverse position, color and polarization. Barred quantities here and below indicate the quantities that are rotated by the R matrix. For a three-gluon Fock state of the three gluon jump operator to the evolution of the density matrix iŝ (4.24) This expression has in principle nine terms. However, for terms that involve the same two gluons connected to a two gluon emission vertex both inM 3 andM † 3 , integration over rapidity produces higher than linear powers in ∆y. For example This term therefore does not contribute to the differential form of the evolution. On the other hand, for the two-gluon-emission vertexes connected to different pairs of gluons, only one explicit ∆y factor arises (4.26) These terms do contribute.
The contribution of the three gluon jump operators is thuŝ (4.27) In the last line we have used superscripts L and R to indicate the position of various factors relative to the density matrixρ v , thusΛ 0,L indicates that this factorΛ 0 is placed to the left ofρ v etc. The ordering is important as the various operators do not commute withρ v . The reason to write the expression in this particular way is that we can use convenient matrix notations, so that products in eq. (4.27) are matrix products over all indexes carried byΛ 0 andb, i.e. color, polarization and transverse coordinate. This pattern clearly generalizes to any odd n. Linear in ∆y contributions arise only from terms where no two gluons are emitted from the same two gluon emission vertex both in M n and M † n . Diagrammatically the terms that yield linear in ∆y contributions are depicted in figure 1.
Generalizing the above analysis to jump operators with 2m + 1 numbers of gluons we obtain  Here i 1 , i 2 , . . . , i 2m+1 is a permutation of 1, 2, . . . , 2m + 1. The sum over P goes over all the possible permutations. The action ofM 2m+1 on the density matrix after summing over all the possible Fock states with 2m + 1 gluons and performing the rapidity integrations, becomeŝ

JHEP05(2020)036
Now adding all the jump operators associated with odd numbers of gluons, their action on the density matrix is

Jump operators with even number of gluons
If the number of gluons in the Fock state is even, the gluons can either be emitted from a two-gluon-emission vertex or from even number of single-gluon-emission vertexes. Since gluons emitted from single-gluon-emission vertexes are uncorrelated in rapidity, the more single gluon emission vertexes are involved, the higher power of ∆y is generated. Recall that we only need to keep terms linear in ∆y. To extract these terms, we allow either no gluons or two gluons to be emitted from the single-gluon-emission vertexes. The rest of the contributions, as we wil see are subleading in powers of ∆y.
The expression for the jump operators with even numbers of gluons byM 2n follows from eq. (4.15) Just like in the case of odd number of gluons, not all the terms in eq. (4.32) contribute to differential evolution. The last term in eq. (4.32) contains an overall factor of (∆y) 2 and therefore can be discarded. The first term in eq. (4.32) does contain terms that are only linear in ∆y, however in the dense limit it is suppressed by a power of α s relative to the second and third terms, as in the dense limit b ∼ 1/g, while Λ ∼ 1. Similar terms arise in the expansion of the normalization factor N (Λ), which we have neglected above. We therefore discard these terms in the week coupling limit. Only the second and third terms in eq. (4.32) are to be evaluated and contribute to the differential evolution of the density matrix. We first evaluateM 0 2n . For example, for a two-gluon Fock state |2 {α 1 ,w 1 ,η 1 ;α 2 ,w 2 ,

JHEP05(2020)036
Generalization to Fock states with 2m gluons is straightforward (4.34) The summation is over all the permutations. The prefactor (1/2) m is needed to account for the fact that Λ is a symmetric matrix. The factor 1/m! takes care of the fact that two permutations that differ only by ordering of some pairs of indices and nothing else, give identical contributions which only need to be counted once.
To calculateM 2 2n we start with the simple situation of two-gluon Fock statê Generalization to Fock states with 2m gluons leads tô (4.36) The action of the two gluon jump operator on the density matrix iŝ For the four gluon operator we similarly find Summing up all the jump operators with even numbers of gluons, we get and its complex conjugate

All together now
We now put together the above results. The wave function renormalization operator and the jump operators with even numbers of gluons contribute to what one might call "the virtual part" of the evolution: (4.41) Adding contributions from jump operators with odd numbers of gluons we get

Operator ordering
As we have mentioned earlier, the relative ordering of operators of b and Λ 0 in eq. (4.42) is not important. It is however important to keep track of the ordering of the classical fields b L and b R relative to the density matrix. More precisely, in the JIMWLK limit one can change the order of various factors in eq. (4.42) as long as each factor (b L −b R ) is kept as a unit and is commuted with any other operator in question. The argument for that was given in [13], and we reproduce it here for completeness.
Recall that the JIMWLK limit is obtained when parametrically b ∼ O(1/g) and Λ 0 ∼ O(1), and additionally the evolution equation should be expanded to order α s . The latter expansion gives the leading contribution when the dense hadron scatters on a dilute target.
Since both b and Λ 0 are functions of j, the commutator between b and Λ 0 can be estimated as The differenceb L −b R can also be estimated as and barred quantities are color rotated byR ∼ O(1).The difference between Λ 0,L and Λ 0,R is estimated similarly

JHEP05(2020)036
Since the factor (b L −b R ) 2 on the right hand side of eq. (4.42) is already of order of α s , the ordering between b and Λ 0 and the difference between Λ 0,L and Λ 0,R contribute to higher orders in α s and therefore can be ignored in the JIMWLK limit as long as one does not order the operators differently in the terms containing b L and b R . Thus for example, one can substitute in eq. (4.42) as a whole, without breaking the factor (b L −b R ) 2 into separate pieces.

The Lindblad form, finally
We now simplify eq. (4.42) using the results of [13]. First we note that the function of Λ appearing in eq. (4.42) can be represented as a square if we indeed forget about the difference between Λ L and Λ R . Define formally We can then write In fact the matrix Θ appears naturally in the calculation of [13]. Since the soft gluon vacuum is a squeezed state due to the presence of the Bogoliubov operatorB, it is the Fock space vacuum of the Bogoliubov transformed set of creation and annihilation operators, related to the original gluon operators a † and a bŷ where Θ and Φ are constrained by the unitarity condition ΘΘ † − ΦΦ † = 1 . (4.52) The following relation was also derived in [13] Using these relations we find

JHEP05(2020)036
After integration over rapidities, eq. (4.54) becomes where N ⊥ = dη(Θ − Φ) has been calculated in [13] where the covariant derivative is defined as D ab i = δ ab ∂ i − igT e ba b e i . Substituting eq. (4.56) into eq. (4.42) , we obtain To write this explicitly as an operator equation we need to make the choice of whether to place the factorN †N to the right or to the left of the density matrix, as both choices, as well as some others are equivalent in the JIMWLK limit. Additionally we need to specify the ordering between the operators b and Λ on one hand and factors of R on another, which has been scrambled in the calculation above. It is not our goal here to carefully restore the correct ordering, as only the JIMWLK limit of this expression is strictly speaking under control. We therefore simply choose the specific ordering which reproduces JIMWLK as well as gives an evolution equation which preserves the Hermiticity of the density matrix also away from the JIMWLK limit. We takeN † where the operator ordering is now specified. With these definitions we write the evolution as Written out explicitly Several words on the nature of the evolution ofρ v as given by (4.61). As we discussed earlier, the only relevant characteristic of a state in the valence Hilbert space is its representation of the color SU(N ) (at each spatial point). Therefore the valence Hilbert space on whichρ v is defined is a direct sum of all the possible subspaces labelled by different representations of SU(N ), i.e. the values of all the Casimir operators at each point JHEP05(2020)036 . . .}, so that H = ⊕H J . Since the density matrix itself depends only on j a , it is a block diagonal operator on this Hilbert space and has nonvanishing matrix elements only between states that belong to the same representation J . The same is true for the operators b α and N αβ since they also are functions of j a only. Thus if not for the operator R in eq. (4.61) the evolution would mix the matrix elements ofρ v in a given representation J only between themselves. The presence of the operatorR changes the nature of the evolution. It shifts j a to j a + gT a and therefore mixes matrix elements ofρ v in one representation with those in another representation with an additional adjoint added in. The operatorR is thus the only source of communication between subspaces of different J in the evolution. Note that even though such cross talk between different representations exists throughout the evolution, the matrixρ v remains block diagonal if it was chosen to be block diagonal at the initial rapidity, since for such an initial condition the right hand side of eq. (4.61) is a function of j a .
Note that eq. (4.61) is somewhat more general than the JIMWLK equation. To obtain the original JIMWLK equation (apart from invoking quantum-classical correspondence which is the subject of the next section) one has to expand (4.61) to second order in Φ a . In fact eq. (4.61) as written here contains both, the JIMWLK limit when expanded in Φ a as well as the KLWMIJ limit when expanded in j a . It can therefore be viewed as an interpolating form of the evolution equation for density matrix between the dense and dilute regimes, just like the corresponding equation for the (quasi)probability density functional W [j] in [13].
The expansion to leading order inΦ a can be performed directly in eq. (4.61). One has to be careful however, since apart from expanding the explicit dependence onΦ a in operators R one also needs to expand the commutators ofρ and b. This is easier done in a somewhat roundabout way, namely transforming the operator equation into the equation for the quasi probability function, performing the expansion there, and then returning to the operator equation using the Wigner-Weyl transformation. We will do precisely this in the next section after introducing the quasiclassical correspondence between the quantum dynamics and dynamics on classical phase space. Here we only present the result of this exercise. The evolution equation in the JIMWLK limit turns out to be and the matrix U is defined as with the path C starting from infinity on the transverse plane and ending at some point x ⊥ .

JHEP05(2020)036
We note that this is precisely the equation for density matrix proposed in [9]. A comment is in order on the form of the evolution equation. First we note that eq. (4.62) is of the Lindblad type. Thus we find that both in the dilute and dense limits the CGC density matrix satisfies Lindblad type equations. Interestingly however, our interpolating equation eq. (4.60), (4.61) does not have the Lindblad form. One might wonder if this absence of Lindblad form is simply an artifact of our approximation. After all, the calculations that lead to eq. (4.60) are under full control only in the two limits. However the reason for deviation from Lindblad in rapidity evolution seems to be very general, and it rather looks like very special conditions have to be satisfied in order for Lindblad form to hold.
We drew an analogy from open quantum systems by treating the valence gluons as "the system" and the soft guons as "the environment". However, time evolution and rapidity evolution are very different concepts. In the former situation, the system degrees of freedom and the environment degrees of freedom are well specified from the beginning and do not change over time. The interaction between the system and the environment is assumed to be Markovian which holds if the system degrees of freedom are slow while the environment degrees of freedom are fast. Time correlations of the environment degrees of freedom are assumed to be local in time compared to the long time scale on which the changes of the system occur, and this leads to Lindblad form of the evolution equation via eq. (2.4).
In the case of rapidity evolution, however, the separation between the "environment"the fast soft gluonic degrees of freedom and the "system" -the slow valence partons is not fixed, but instead the separation boundary moves together with the evolution parameter. As the rapidity increases one therefore integrates over additional degrees of freedom, namely those whose rapidity label is between the old and new values of the evolution parameter -∆y. The rapidity thus appears not only as a parameter of the evolution analogous to time, but also as the label of the quantum states which are being integrated out in the process. This integration over additional degrees of freedom is part and parcel of rapidity evolution, and the increment in the density matrix ∆ρ v proportional to ∆y arises due to this integration. Thus the "Markovian" regime (i.e. short correlation length of the environment modes in rapidity) albeit sufficient to guarantee existence of differential evolution, does not guarantee that this evolution is of Lindblad form. In fact it is easy to trace that the additional integration over the rapidity label is in fact the reason why the general argument that leads to Lindblad form for time evolution in quantum mechanics, is violated in our calculation.
The crucial missing piece is eq. (2.4). As discussed in section II, in an evolving quantum system for small ∆t the probabilityM † n (dt)M n (dt) ∝ dt is proportional to dt for every environment state n save the vacuum. One can then define the jump operatorL n viâ M n (dt) = √ dtL n and the Lindblad form follows. On the other hand, in the case of rapidity evolution the factor ∆y arises only as a result of the integration over rapidity label of the gluon states in the rapidity window [y ; y + ∆y]. As a result we have but the probabilities for individual states with fixed η do not scale with ∆y. It thus does not appear to be possible in general to define a jump operator unlessM n (η) has very special properties. Indeed, examining our calculation, for example in eq. (4.32), we realize that the fact that only the termsM 0 2nρ vM 2 † 2n +M 2 2nρ vM 0 † 2n contribute at linear order in ∆y is precisely due to the integration over the rapidities of the 2n gluons. This is also the reason this contribution cannot be written in the standard Lindblad form ∆yL 2nρvL † 2n . The same is evidently true also for the odd gluon contributions.
Nevertheless in some special cases the Lindblad form may be attainable. For example, ifM (η, ∆y) =M (∆y), i.e. if the probability of a particular state does not depend on the rapidity of the gluon, the jump operator can indeed be defined. This is precisely the situation we encounter in the derivation in the dilute (KLWMIJ) limit. In this case the coherent operator C involves only the gluon creation operator integrated over rapidity and as a result the probabilityM † n M n does not depend on the rapidity label of the gluons in the state n. This then allows to take the "square root" of the probability and define the corresponding jump operatorL n which ensures that evolution is in Lindblad form. It is more difficult to trace the origin of the Lindblad form in the dense limit. However, given that JIMWLK and KLWMIJ limits are dual to each other, it is not surprising that such a form indeed exists.
In this section we have discussed the energy evolution in terms of the CGC density matrix. This is not the way it has been formulated in the literature so far. In the next section we show how to relate the two formulations.

From the Lindblad equation to the JIMWLK equation via quantumclassical correspondence
In the previous section we have derived the rapidity evolution equation for the CGC density matrix. To turn this into the conventional JIMWLK evolution equation we will invoke a variant of the Quantum-Classical correspondence, which for simple quantum mechanical systems has been studied 50 years ago, for review see [14]. To start with, we present this analysis as it appears in [14,[24][25][26].

Quantum mechanics in phase space
For simplicity, consider a system with one degree of freedom equipped with the canonical variablesq andp that satisfy the cannonical commutation realtion [q,p] = i. This can be achieved if one can find a one to one correspondence between an arbitrary quantum operator A(q,p) and a corresponding classical function A(q, p), and additionally similar correspondence for the density matrixρ → W (q, p).

JHEP05(2020)036
In principle one can devise different mappings that achieve this goal. One widely used mapping is the Wigner-Weyl transformation that maps fully symmetrized quantum operators in the Hilbert space to the corresponding classical functions in the phase space (and back). The familiar form of the Wigner transformation uses eigenstates of theq operator and defines the classical phase space functions via The Wigner function W [q, p] that corresponds classically to the density matrix is often called the quasi probability distribution function on the phase space. To formulate the mapping in a basis-independent way, one follows Weyl's correspondence rule which associates fully symmetrized operators in Hilbert space to classical functions in phase space. Consider the following representation of an operator G(p,q) This can be regarded as an operator Fourier transformation. First of all, note that the operator G written in this form is necessarily symmetric under permutations ofq andp. This is straightforward to see by expanding the exponential in Taylor series. This is however not a restriction on the set of operators one can consider, as any operator function can be written in a symmetric form utilizing the commutation relation betweenp andq. The simplest example of such symmetrization ispq = 1 2 (pq +qp) − i 2 . One can easily convince oneself that any polynomial ofp andq can be written in a symmetric form of this type.
Given the representation eq. (5.4) we define a classical function on the phase space via This is Weyl's rule for correspondence between quantum operators and classical functions on phase space. Note that under Weyl's rule, the same Fourier kernel f (u, v) is used in eqs (5.5) and (5.4). From eqs. (5.5) and (5.4), the mapping between F (p, q) and G[p,q] can be represented as

JHEP05(2020)036
The Weyl mapping kernel ∆ w (p −p, q −q) is a functional of both canonical operatorsp,q and phase space classical variables q, p. It has the following properties The basis independent mapping in eqs. (5.6) reproduces the familiar Wigner transformation eq. (5.3) when the eigenbasis ofq is used.
where we have used the Baker-Campbell-Hausdorff formula, and Thus the Weyl's correspondence rule provides a basis independent mapping between quantum operators in Hilbert space and classical functions in phase space. Using eq. (5.6) one finds that expectation values of quantum observables can be calculated as weighted integrals in the phase space with W (p, q) = Tr ρ∆ w (p −p, q −q) . (5.13) Note that once the operator A s [q,p] is written in a fully symmetrized form with respect toq,p, its associated Wigner-Weyl mapped classical function can be obtained by simply replacing the quantum operatorsq,p with their classical counterparts q, p, Consider now a product of two operatorsF [q,p] =Â[q,p]B[q,p]. Even if bothÂ[q,p] andB[q,p] are fully symmetrized, their product as a function ofp andq does not necessarily have a fully symmetrized form, and therefore the Wigner-Weyl transform of a product is not equal to product of two Wigner-Weyl transforms. Instead, the correct procedure to obtain the Wigner-Weyl transformation for a product of two observables is where the derivatives act on functions on the left or on the right as indicated by the arrows.

JHEP05(2020)036
An alternative representation for the transformation of products of operators can be achieved by introducing the left and right Bopp operators and The two sets of Bopp operators labelled by "L" and "R" are operators in phase space rather than Hilbert space, as they act on classical functions of p and q. Note that the pairs (Q L , P L ) and (Q R , −P R ) as operators in phase space form the same Heisenberg algebra as do (q,p) in the Hilbert space. Using Bopp operators, the Wigner-Weyl transformation of a product of two observables is expressed as Here A s [Q L , P L ] is obtained by replacingq andp in A s [q,p] by Q L and P L , respectively, and similarly for B[Q R , P R ]. One can use either set of Bopp operators, depending on whether one uses them in the left factor or the right factor of the product. The two expressions in eq. (5.17) are equivalent. The action of Bopp operators represents the additional symmetrization rearrangement necessary in order to represent a product of two symmetrized operators in a completely symmetrized form.
As an application of the

Quantum-classical correspondence for SU(N ) charges
In the context of the JIMWLK evolution, the relation between the "probability density functional" W [j] and the density matrixρ is similar to that between W [p, q] andρ in quantum mechanics as described in the previous subsection. In this subsection we make this relation explicit. Much of this discussion already appears in the literature, e.g. in [10], but the relation to the classical-quantum correspondence and Wigner-Weyl transformation has not been elucidated in the past. Unlike the canonical case discussed above, we are now dealing with the system whose phase space is spanned by the generators of the SU(N ) group j a (x ⊥ ). It is important to stress that the components of color charge density are coordinates on the phase space, and not on configurations space. The Hilbert space of the corresponding quantum system is spanned by the quantum operatorsĵ a (x ⊥ ) that satisfy the SU(N ) commutation relations [ĵ a (x ⊥ ,ĵ b (y ⊥ )] = igf abcĵc (x ⊥ )δ(x ⊥ − y ⊥ ). Note that although the full Hilbert space of CGC requires introduction of the operatorsΦ a (x ⊥ ), the observables that are currently considered in all calculations are only functions of the color charge densityĵ a (x ⊥ ). It is thus sufficient for our purposes to discuss the quantum-classical correspondence for operators that depend only on j(x ⊥ ). One must keep in mind however that if one wishes to generalize the framework along the lines of [9], this correspondence has to be extended to include also functions ofΦ(x ⊥ ).
For quantum systems of spins, most notably the SU(2) group, the mapping between the Hilbert space and the phase space has long been studied [27][28][29][30]. These studies mostly rely on introducing a particular (over)complete basis (ususally generalized coherent states) and working in a fixed representation of the underlying Lie group. Our situation is slightly different, since as discussed above the valence Hilbert space is a direct sum of different representations of SU(N ). We therefore cannot fix the representation and instead will rely on the operator properties of the quantum-classical correspondence. The purpose of this section is, drawing analogy to the canonical case to provide the mapping between quantum operators in valence Hilbert space and classical variables in non-Abelian phase space, as well as relation between the quantum density matrix and "classical" quasi probability distribution.
We concentrate on operators in the Hilbert space which can be written as functions ofĵ a (x) and that are fully symmetric with respect to interchange of the different color components ofĵ a . If an operator is not written in a fully symmetrized form, it can always be brought into such form by repeated use of the basic commutation relations ofĵ a . This has been explicitly demonstrated in [10].
We will construct the analog of Weyl's quantum-classical correspondence where the classical counterpart of a quantum operator is obtained by replacing an operatorĵ a with its classical counterpart j a once a quantum operatorÔ is expressed in terms of fully symmetrized products ofĵ a i , so thatÔ = O s (ĵ a ), where the subscript "s" indicates a fully symmetric function. This is the most straightforward generalization of Weyl's quantumclassical correspondence. In the following the spatial coordinates ofĵ a (x ⊥ ) are suppressed JHEP05(2020)036 sinceĵ a (x ⊥ ) with different transverse coordinates commute. All the nontrivial action therefore happens at the same transverse coordinate.
In analogy with the discussion of the Weyl's correspondence rules for canonical operators in the previous subsection, we adopt the following rules for correspondence between an operator and a classical function on the phase space (5.23) Note that for SU(N c ), the color index runs from 1 to N 2 c − 1. The above Fourier transformations are understood as N 2 c − 1-variate transformations. It is easy to see by Taylor expanding the second of eq. (5.23) that G s is a fully symmetric function ofĵ a .
Note that eq. (5.23) is an operator relation, and is not limited to any particular representation of the SU(N ) group, but is rather valid on all the valence Hilbert space.
Eq. (5.23) leads to the following relation between G s [ĵ] and F [j] with the mapping kernel Our definition is such that the integration dj is over all real valued j a with a simple integration measure on R N 2 −1 . The ∆ W (j,ĵ) maps classical functions F [j] to quantum operators G s [ĵ a ]. By requiring that for a fully symmetrized operator G s [ĵ] the corresponding classical function in phase space is just G s [j], so that F [j] = G s [j] we can also find the inverse mapping.
Let us write this mapping in the suggestive form: Substituting eq. (5.24) into eq. (5.26), one obtains the condition that∆ W [j,ĵ] must satisfy The following expression of∆ W [j,ĵ] solves this constraint: Here dg α is the Haar measure over the SU(N ) group. Each group element g α is labelled by the parameters α. The factor d r (ĵ) denotes the dimension of the particular representation r and is viewed here as a function ofĵ which depends only on the Casimir operators of the Lie algebra. The function is such that for a given representation its numerical value is equal to the dimension of this representation.

JHEP05(2020)036
To prove that∆ W [j,ĵ] satisfies eq. (5.26), we use the Peter-Weyl theorem [31] for representations of Lie group ∞ r=1 dr j,k=1 . (5.29) Here D (r) (g α ) is the representation of group element g α and r indicates all the irreducible, inequivalent representations. As above, d r is the dimension of the representation r. Using eq. (5.28) in the left hand side of eq. (5.27), and remembering that summation over all representations r is a part of tracing over the valence Hilbert space, we recover the right hand side of eq. (5.27).
The above expression is the formal definition of the kernel∆, however for all practical purposes one does not need to know its explicit form. This is because the classical counterpart of the operator G[ĵ] is simply obtained by substitutionĵ → j once the operatorĜ is written in the fully symmetrized form. Now we can establish the relation between quantum average of operators in Hilbert space and phase space weighted integrations where the summation goes over all possible permutations of (1, . . . , n).
The other issue we need to understand in order to formulate evolution in the classical phase space approach is how to extend the mapping for products of quantum operators. In principle, this involves generalizing Moyal's star-product to a general Lie algebra. However, rather than taking this general mathematical approach, the particular realization for SU(N ) algebra has been worked out in [10,32]. Consider a product of two operators with each one written in the symmetrized form  The Bopp operators j a L and j b R act on functions in the phase space rather than the Hilbert space. It is straightforward if somewhat tedious to explicitly check that, similarly to Bopp operators defined in eqs. (5.15), (5.16), the SU(N ) phase space Bopp operators j a L and −j a R form the same SU(N ) algebra as the operatorsĵ a on the Hilbert space. In addition, j a L and j b R commute [j a L , j b R ] = 0. We find it interesting to note that the functional form of the Bopp operators j a R involves exactly the same function as in eq. (3.19), which ensured correct operator properties of the charge shift operatorR.
As a corollary to this discussion consider Hermitian conjugation in Hilbert space x ⊥ ], depending on its position relative to the factorρ in eq. (4.61). The operator N †

JHEP05(2020)036
Additionally we need to understand how the operatorR is mapped to the phase space. To do this, we note that for a fully symmetrized operatorÔ s [ĵ a ], the action ofR operator iŝ with the phase space shift operator Therefore, we should simply replace the action of the operatorR with the phase space shift operator R p . It is now straightforward to write the equation for W . We obtain One can now take various limits to reproduce the results known in the literature. In particular, assuming that j is small and expanding the right hand side of eq. (5.42) to second order in j one straightforwardly recovers the so called KLWMIJ equation [23,33].
Alternatively, keeping all orders in j, but expanding to second order in δ/δj one reproduces the JIMWLK equation. This last expansion is a little more involved, but it is performed explicitly in [22]. We reproduce the derivation here for completeness.
To reproduce the JIMWLK kernel, we truncate R p to first order in δ/δj a and expand b L and b R around b(j a ). We only need to keep first order terms, since eq. (5.42) contains a factor (b L − b R ) 2 .
At this order there is no need to expand N ⊥,R or N ⊥,L so that both are substituted by N ⊥ (j). We have .
In the last line we have used igT e ba b e i = δ ab ∂ i − D ab i and igT e ba j e = −(∂D − D∂) ab = igj e T a eb as well as be (x ⊥ , z ⊥ ). Additionally using eq. (4.57), we calculate , (5.44) where the standard eikonal scattering matrix U is defined in eq. (4.64).

JHEP05(2020)036
Eq. (5.42) now becomes One can express the JIMWLK equation in terms of single gluon scattering matrix U (x ⊥ ) ab rather than the color charge density j a (x ⊥ ). The relation between the two was derived in [22]. Now we can relate functional derivatives with respect to the U matrix to those with respect to the color current. (5.47) In the above we have used U † T a U = U ab T b and D l = U † ∂ l U . We have also used
(5.50) Also recall that The standard JIMWLK kernel is reproduced as (5.53)

From JIMWLK back to Lindblad
Finally using eq. (5.45) and eq. (5.44) we can transform the evolution equation back to Hilbert space. First we note, that the amplitude Q a i is hermitian as an operator on the phase space. This is obvious from eq. (5.52), since J b R (y ⊥ ) can be commuted to the left through the factors of U , as the right color index on U is contracted with the index of J b R (y ⊥ ). One can then write the JIMWLK equation as where the commutator is understood as the commutator of the operators on phase space. Transforming this back to Hilbert space we see that to this order all we need to do is substitute −i δ δj a (z) →Φ a (z) and keep the structure of the double commutator. This procedure gives eq. (4.62) as claimed.

Discussion
This paper is devoted to analysis of the high energy limit of hadronic scattering, and its energy evolution formulated as effective quantum theory.
The dynamics of this effective theory is governed by a density matrix. Here we were able to define this "reduced" density matrix in a way reminiscent to an open quantum system, i.e. bi partitioning the degrees of freedom into the "system" and "environment" and integrating over the environment. In the present case the bi partitioning is into the "valence" gluons as the system and "soft" gluons as the environment. Despite some similarities, there is a significant difference between the high energy limit considered here and bi partitioning in a quantum open system. In a quantum system one normally integrates JHEP05(2020)036 completely the environment and then considers only observables that depend on the degrees of freedom of the "system". This is not the case for the high energy limit, since the soft gluons contribute nontrivially to the color charge density, which is the basic observable in the effective theory. Defining the reduced density matrix is therefore rather nontrivial. Nevertheless we were able to do it.
We have then followed the usual assumption made in the derivation of the high energy evolution, i.e. that only the distribution of the color charge density in the transverse plane is relevant for determining hadronic properties at high energy. Assuming that the reduced density matrix of a hadronic system depends only on the components of the color charge density results in a quasi diagonal density matrix in the sense that its matrix elements between states belonging to different color representation vanish. This is true both in dense and dilute limits, and generalizes the notion of diagonal density matrix discussed in [9] to arbitrary parametric values of color charge density.
Under this assumption we have shown that the rapidity evolution of the reduced density matrix is of the Lindblad type in the two limiting cases -the dense (JIMWLK) and the dilute (KLWMIJ) limits. This is true even though the nature of the energy evolution is in principle quite different from the nature of time evolution of a dynamical quantum system.
Interestingly the evolution equation that interpolates between the two limits, eq. (4.60) does not have a Lindblad form. Although the derivation of this interpolating equation is not under parametric control, the basic features of the derivation are generic, and our analysis shows that the absence of Lindblad form should be a rule rather than exception. The basic reason is that the rapidity plays a dual role in high energy evolution: it is the analog of the evolution time on one hand, and is a quantum number that labels the quantum states of the environment that are integrated out on the other hand. This invalidates in principle the usual argument for the Lindblad form of the differential evolution equation.
Finally, we have shown how to rigorously relate the reduced density matrix description of the evolution with the approach used in most pertinent literature based on the probability density functional W [j]. To this end we have explored the Wigner-Weyl transformation, which maps the Hilbert space description of a quantum system in terms of density matrix ρ into the classical phase space description in terms of quasi probability distribution W . By adapting this transformation to the present case we have shown that the Lindblad evolution equation forρ is indeed equivalent to a Fokker-Planck type equation for W , where the components of the color charge density j are considered as coordinates on a classical non Abelian phase space. This Fokker-Planck equation reduces to KLWMIJ and JIMWLK equations in the appropriate dilute and dense limits.
In quantum optics, it has been known for decades that the Lindblad master equation maps to the Fokker-Planck equation through quantum-classical correspondence. Here we have established the same in the context of high energy evolution with the JIMWLK (or KLWMIJ) playing the role of the Focker-Planck equation.
We stress again that this paper deals only with the conventional JIMWLK/KLWMIJ setup, where the density matrix is assumed to depend only on the color charge density degrees of freedom. Recently it was suggested that this framework may be too restrictive and may not be adequate for studying some interesting observables at high energy [9]. Such observables, like correlations between the transverse momentum and density in the JHEP05(2020)036 transverse plane may be formally subleading at high energy, but could be of great interest in the study of correlations in particle production. In order to include these into consideration one has to extend the conventional framework and allow for density matrices that depend not just on color charge densityĵ, but also on their conjugate variables, which in the present paper we have identified as the operatorsΦ a . It was suggested in [9] that the evolution of this more general density matrix is also given by the same Lindblad equation. Although this seems very likely to be the case, our current derivation does not cover this interesting more general situation. It should be possible to extend our current method to deal with this intriguing problem. This investigation is currently under way.
We now comment on several questions/issues that arise from our results. First, the fact that the evolution of the reduced density matrix beyond dense-dilute limit is most likely not of Lindblad type begs an interesting general question. It is known that Lindblad equation preserves the properties of the density matrix, namely normalization and positivity. Is this also the case for eq. (4.60) even though it is not in Lindblad form? It is quite obvious that the normalization of the density matrix is preserved under eq. (4.60), since its right hand side is a commutator, and therefore has a vanishing trace. As for the positivity, it is more difficult to establish. We note however, that the differential evolution follows from the Krauss representation eq. (3.23) which does preserve positivity [15]. We therefore believe that the differential evolution eq. (4.60) does indeed preserve positivity and thus is a consistent evolution of a density matrix. If this is the case, one is lead to a general conclusion that the set of possible differential evolutions of a density matrix is not limited to equations of Lindblad type.
Second, we note that one of the useful perspectives on the JIMWLK evolution is that of a Langevin equation for Brownian motion in the space of Wilson line U ab (x ⊥ ) [34,35]. The bi partitioning into the "system" of the hard gluons and "environment" due to the soft gluons harmonizes nicely with the random walk picture. After a boost by ∆y, the soft gluons can be emitted into any of the multigluon Fock states {|n }. This emission contributes to a random addition to the color charge density j a y+∆y ∼ j a y + δj a with δj a being a random variable, which therefore random walks in the color space. The Langevin equation is a reformulation of the Focker-Planck equation, which is equivalent to JIMWLK. It is then interesting to ask whether such a Langevin description can be extended beyond the leading order. The NLO JIMWLK equation has been derived some years ago [36][37][38][39]. Naturally the derivation involves integration over gluon and quark degrees of freedom in the rapidity interval ∆y. As opposed to the leading order, where as we discussed the probability to create all single gluon states is equal, independent of their rapidities, at NLO there is a genuine integration over rapidity of two soft parton states. This suggests that the evolution equation for the density matrix is not of Lindblad type for the same reason eq. (4.60) is not. If that is the case one does not expect it to be equivalent to a Focker-Planck equation for the quasi probability function and thus the Langevin description may well not be possible.
We hope that the new perspective on high energy evolution discussed in this paper will be useful not only for a more fundamental understanding of JIMWLK equation but will also prove useful for future developments.

JHEP05(2020)036
A The operatorΦ a In this appendix we present the derivation of the quantum "phase" operatorΦ a , which via eq. (3.8) defines the quantum shift operator of the color charge density. Part of this derivation appears in the text, but we keep here all the details for completeness.
We are looking for M ab (Φ), as a functinal of operatorΦ a that satisfies the commutation relations  Here we have used −if abc = T a bc . From the above explicit calculations, it is natural to make the following ansatz We here check the consistency condition eq. (A.15) order by order to verify that it indeed is satisfied. We do not have a proof for the case of general N, but we believe the same procedure can be carried out to high order terms. The left hand side of eq. (A.15) has N = m + n − 1 in terms of power ofΦ a . Also note that m ≥ 0 and n ≥ 1. In the following, we calculate the cases N = 0, 1, 2, 3 in details.
The guiding principles of organizing these terms are the following: • d 1 , d 2 , d 3 are symmetric indices, we are free to interchange among them.
• We rearrange terms according to their b, c indices. If b, c are indices of the same matrix like T m bc , there is no need for further simplication. If we have T b , T c in adjacent position like T b T c , the b ↔ c subtraction gives a commutator, which results in placing b, c indices on the same matrix T . Then no further simplification is needed.
• If T b , T c are not directly in the adjacent position, we rearrange the b, c as the indices of matrix element not the label of T matrix.
• We use the cancellations between terms like (T a T d 1 T d 2 T d 3 ) bc and (T d 1 T d 2 T d 3 T a ) cb .

JHEP05(2020)036
First note the last two terms in the second bracket .25) and (A. 26) the first two terms in the second bracket In obtaining the third equality, we have moved the T a matrix in the first term gradually to the far right so that the resulting term will cancel the third term. We also moved the T a matrix in the second term passing T d 2 , which then will cancel the fourth term. Now consider the terms in the first bracket.
Note that c 2 2 = (1/12) 2 = 5 6! and c 0 c 4 = − 1 6! . Therefore for N = 3 terms eq. (A.15) holds. We have checked four leading terms in the expansion of the Jacobi identity. The same explicit procedure can be followed for higher order terms as well, but it becomes increasingly cumbersome. We therefore stop at this point.