Next-to-leading order Balitsky-Kovchegov equation beyond large Nc

This is a self-archived version of an original article. This version may differ from the original in pagination and typographic details. Author(s): Title: Year: Version: Copyright: Rights: Rights url: Please cite the original version: CC BY 4.0 https://creativecommons.org/licenses/by/4.0/ Next-to-leading order Balitsky-Kovchegov equation beyond large Nc © Authors, 2020 Published version Lappi, Tuomas; Mäntysaari, Heikki; Ramnath, Andrecia Lappi, T., Mäntysaari, H., & Ramnath, A. (2020). Next-to-leading order Balitsky-Kovchegov equation beyond large Nc. Physical Review D, 102(7), Article 074027. https://doi.org/10.1103/physrevd.102.074027 2020 Next-to-leading order Balitsky-Kovchegov equation beyond large Nc T. Lappi , H. Mäntysaari , and A. Ramnath Department of Physics, University of Jyväskylä, P.O. Box 35, 40014 University of Jyväskylä, Finland and Helsinki Institute of Physics, P.O. Box 64, 00014 University of Helsinki, Finland (Received 12 July 2020; accepted 2 October 2020; published 29 October 2020) We calculate finite-Nccorrections to the next-to-leading order (NLO) Balitsky-Kovchegov (BK) equation. We find analytical expressions for the necessary correlators of six Wilson lines in terms of the two-point function using the Gaussian approximation. In a suitable basis, the problem reduces from the diagonalization of a six-by-six matrix to the diagonalization of a three-by-three matrix, which can easily be done analytically. We study numerically the effects of these finite-Nccorrections on the NLO BK equation. In general, we find that the finite-Nccorrections are smaller than the expected 1=N2 c∼10%. The corrections may be large for individual correlators, but have less of an influence on the shape of the amplitude as a function of the dipole size. They have an even smaller effect on the evolution speed as a function of rapidity. DOI: 10.1103/PhysRevD.102.074027 I. INTRODUCTION In hadronic collisions at high energies, large gluon densities are created by the emission of soft gluons carrying a smallfractionofthelongitudinal momentumoftheparent[1]. Nonlinear dynamics of gluons becomes important in such an environment, where parton densities eventually grow to become on the order of the inverse of the QCD coupling αs. To describe QCD in thisregion, thecolor glass condensate (CGC) effective field theory [2] has been developed. In the CGC framework, cross sections for various scattering processes can be expressed in terms of correlators of Wilson lines. A Wilson line describes the eikonal propagation of a parton in the strong color field of the target. The energy dependence of the target color fields, and thus cross sections, is obtained by solving the socalled Jalilian-Marian–Iancu–McLerran–Weigert–Leonidov– Kovner (JIMWLK) equation [3–6].Thisisaperturbative evolution equation that describes the Bjorken-xdependence of a Wilson line. In phenomenological applications, it is usually convenient to work directly in terms of the Wilson line correlators, and to solve instead the Balitsky-Kovchegov (BK) equation [7,8] for the dipole operator (correlator of two Wilson lines), which can be obtained from the JIMWLK equation in the large-Nclimit. The CGC framework has been used extensively in phenomenological applications at leading order (LO) in αs, with the evolution equations resumming contributions ∼αsln 1=x to all orders. Running coupling effects derived in Refs. [9–12] (see also [13]) can also be taken into account. The nonperturbative initial condition for the smallxevolution is obtained by performing fits to the HERA structure function data [1,14], for example in Refs. [15–18] (see also [19,20]). The obtained initial condition can then be used for various calculations, for example particle production in proton-nucleus collisions [17,21–27].In the future, the nuclear deep inelastic scattering (DIS) experiments at the Electron Ion Collider (EIC) [28,29] in the US, at the LHeC [30] at CERN and at the EicC in China [31] will provide a vast amount of precise data from clean DIS processes. These experiments will be able to probe the nuclear structure where nonlinearities are enhanced by roughly A1=3higher densities compared to the proton. Before the EIC, similar studies limited to the photoproduction region can be performed in ultra-peripheral heavyion collisions [32,33]. In order to quantitatively study nonlinear dynamics in high-energy scattering processes (and especially at the future EIC), it is crucial to move beyond LO accuracy. The next-to-leading order (NLO) evolution equations are available:the NLOBKequationwas derived inRef. [34]and the NLO JIMWLK equation was derived in Refs. [35,36]. Similarly, the impact factors are becoming available at NLO for some processes: inclusive DIS [37–41] (in the case of massless quarks), exclusive vector meson production [42,43] (see also [44]) and particle production in proton– nucleus collisions [45]. However, the phenomenological applications of these are still developing [46–51]. The BK equation is usually solved in the large-Nclimit. In the LO case, the large-Nclimit makes it possible to Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article’s title, journal citation, and DOI. Funded by SCOAP3. PHYSICAL REVIEW D 102, 074027 (2020) Editors' Suggestion 2470-0010=2020=102(7)=074027(18) 074027-1 Published by the American Physical Society express the four-point correlator of fundamental representation Wilson lines in terms of the two-point function. In detailed numerical studies, it has been shown that the finiteNccorrections are smaller than the naive expectation of Oð1=N2 cÞ[52,53]. At NLO, the equation involves six-point functions of fundamental Wilson lines that must similarly be expressed in terms of the two-point function in order to close the equation. The purpose of this work is to see if the finite-Nccorrections are similarly small in the case of the NLO equation, where all corrections of the order α2 sare taken into account. In order to numerically solve the BK equation at finite Nc, we use the Gaussian approximation [54–56] to derive analytical parametric equations for the six-point correlators in terms of the two-point correlators. We study numerically the finite-Nccorrections to these correlators, and also their effect on the NLO BK evolution. In addition to the BK equation, higher-point correlators are needed in the calculations of multi-particle correlations in the CGC framework, see eg. Refs. [57–60]. The structure of the paper is as follows. In Sec. II, we introduce the NLO BK equation and provide both the large-Ncand finite-Ncexpressions for the correlators that will be studied. In Sec. III, we introduce the Gaussian approximation, explain the diagrammatic notation used in the rest of the paper and then explain the analytical calculation done for finding the parametric equations for the six-point correlators. Section IV contains the numerical results obtained from using the analytical expressions for the six-point correlators to solve the BK equation. Finally, we end with a few concluding remarks and a summary of our main results. II. THE BK EQUATION AT NLO For any product of n=2pairs of fundamental Wilson lines UU†, we use the notation SðnÞ x1;x2;…;xn−1;xn≔1 Nc trðUx1U† x2…Uxn−1U† xnÞ:ð1Þ The NLO BK equation in the case of zero active quark flavors (nf¼0) reads [34] ∂YhSð2Þ x;yi¼αsNc 2π2KBC 1⊗hD1iþα2 sN2 c 16π4K2;1⊗hD2;1iþα2 sN2 c 16π4K2;2⊗hD2;2iþOðnfÞ;ð2Þ where the brackets hi refer to the expectation value over target color field configurations. The kernels are KBC 1¼r2 X2Y21þαsNc 4πβ Nc ln r2μ2−β Nc X2−Y2 r2ln X2 Y2þ67 9−π2 3−10 9 nf Nc −2ln X2 r2ln Y2 r2;ð3Þ K2;1¼−4 Z4þ2X2Y02þX02Y2−4r2Z2 Z4ðX2Y02−X02Y2Þþr4 X2Y02−X02Y21 X2Y02þ1 Y2X02 þr2 Z21 X2Y02−1 X02Y2×lnX2Y02 X02Y2;ð4Þ K2;2¼r2 Z21 X2Y02þ1 Y2X02−r4 X2Y02X02Y2ln X2Y02 X02Y2:ð5Þ The convolutions ⊗in Eq. (2) denote integrations over the transverse coordinate z(in KBC 1)orzand z0(in K2;1and K2;2). We use the notation r2¼ðx−yÞ2,X2¼ðx−zÞ2, X02¼ðx−z0Þ2,Y2¼ðy−zÞ2,Y02¼ðy−z0Þ2and Z2¼ðz−z0Þ2. We note that the kernel proportional to nfis also available [34]. Since the purpose of this work is to study the importance of the finite-Nccorrections in the NLO BK equation, we do not include contributions proportional to nf. The finite-Nceffects could be expected to be similar to the nf¼0case. The Wilson line operators appearing in Eq. (2) are hD1i¼hSð2Þ x;zSð2Þ z;yi−hSð2Þ x;yi;ð6Þ hD2;1i¼hSð2Þ x;zSð2Þ z;z0Sð2Þ z0;yi−1 N2 chSð6Þ x;z;z0;y;z;z0i −ðz0→zÞ;ð7Þ hD2;2i¼hSð2Þ x;zSð2Þ z;z0Sð2Þ z0;yi−ðz0→zÞ:ð8Þ Although the original NLO BK equation in the form presented in Ref. [34] does not contain the subtraction z0→zin D2;2, we have introduced the subtraction to improve numerical stability. This subtraction term has no effect on the final evolution because the integral of K2;2 over z0vanishes if the Wilson line operator term does not T. LAPPI, H. MÄNTYSAARI, and A. RAMNATH PHYS. REV. D 102, 074027 (2020) 074027-2 depend on z0(see Ref. [34]). We will refer throughout this work to two pieces of the right side of Eq. (2) as the (i) αsNc 2π2KBC 1⊗hD1i∼“LO-like”contribution, (ii) α2 sN2 c 16π4K2;1⊗hD2;1iþα2 sN2 c 16π4K2;2⊗hD2;2i∼“NLO-like” contribution. In other words, we separate the terms in the NLO BK equation by the types of Wilson line correlators, not by the order in αs. Thus, the LO-like contribution also includes a significant α2 scorrection. The interpretation of the NLO BK equation is that one considers all possible ways to emit either one or two gluons, at transverse coordinates zand z0, from the boosted dipole consisting of quarks at transverse coordinates xand y. The effect of the boost is that instead of the original dipole projectile, the original quarks and the emitted gluons scatter off the target color field. As such, the evolution can be seen to describe the evolution of the projectile probing the target structure. On the other hand, the emitted gluons can also be taken to be a part of the target wave function, in which case the boost corresponds to the evolving target color field as probed by the original projectile. For a more detailed discussion on the NLO evolution in the projectile or target wave function, the reader is referred to Ref. [61]. The NLO BK equation is known to be unstable [62] due to the large contributions enhanced by the large double transverse logarithm ln X2=r2ln Y2=r2. We resum these contributions to all orders following the procedure developed in Ref. [63], which was numerically confirmed in Ref. [64] to result in a stable evolution (see also Ref. [65] for an equivalent resummation of the same double logarithms). In addition, we include the running of the QCD coupling by noticing that the terms proportional to the beta function coefficient βin Eq. (3) should be resummed into the running coupling. We implement this resummation by following the Balitsky prescription from Ref. [12]. Both running coupling and double transverse logarithm resummations are included by modifying the kernel KBC 1as αsNc 2π2KBC 1→ αsðrÞNc 2π2KDLAr2 X2Y2þ1 X2αsðXÞ αsðYÞ−1 þ1 Y2αsðYÞ αsðXÞ−1þKfin 1:ð9Þ The double log corrections to all orders are taken into account by the factor KDLA ¼J1ð2ffiffiffiffiffiffiffiffiffi ¯ αsx2 pÞ ffiffiffiffiffiffiffiffiffi ¯ αsx2 p;ð10Þ where ¯ αs¼αsNc=π. The double logarithm here is x¼ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ln X2=r2ln Y2=r2 p.IflnX2=r2ln Y2=r2<0, then an absolute value is used and the Bessel function is changed from J1→I1(see Ref. [63]). The scale of the coupling in KDLA is determined by the smallest dipole minfr2;X2;Y2g. In addition to the double log contributions, one can also resum a set of higher-order contributions enhanced by single transverse logarithms, as shown in Ref. [66]. For the purposes of this paper, this resummation is not necessary and is excluded for simplicity. In this running coupling prescription, we keep the order α2 sterms in the kernel KBC 1 that are not proportional to the beta function. These are included in the term Kfin 1, which reads Kfin 1¼α2 sðrÞN2 c 8π3 r2 X2Y267 9−π2 3−10 9 nf Nc:ð11Þ The strong coupling constant in the transverse coordinate space is evaluated as αsðrÞ¼ 4π βlnf½ð μ2 0 Λ2 QCDÞ 1 cþð 4C2 r2Λ2 QCDÞ1 ccg ;ð12Þ where β¼ð11Nc−2nfÞ=Nc. We take nfto be zero in both Kfin 1and β. We use1C2¼1and μ0=ΛQCD ¼2.5in our numerical calculations, which freezes the coupling at αsðr→∞Þ¼0.762 in the infrared, and c¼0.2which controls the transition to the infrared region. The initial condition for the BK equation is taken from the McLerran-Venugopalan (MV) model [67,68]. In the MV model, the color charge density is assumed to be a random Gaussian variable, with a zero expectation value and a variance proportional to the local saturation scale Q2 s. The dipole correlator in the MV model is written as hSð2Þ x;yiMV ¼exp −r2Q2 s0 4ln 1 rΛQCD þe:ð13Þ Here, the constant eacts as an infrared regulator. We use ΛQCD ¼0.241GeV and Q2 s0¼1GeV2in the numerical analysis. In analytical studies of the correlators of Wilson lines in specific “line”coordinate configurations in Sec. IVA,we use the GBW [69] form for the dipole correlator hSð2Þ x;yiGBW ¼exp −r2Q2 s0 4:ð14Þ In principle, the resummation procedure for the double transverse logs would also change the initial condition, as discussed in Refs. [63,64]. However, since the initial condition is a nonperturbative input for the evolution, we consider Eq. (13) to be the nonperturbative initial condition for the resummed evolution as well. For the 1A generic estimate [9,13] would be C2¼e−2γE≈0.32.We use a larger value C2¼1which results in slightly slower evolution, as the C2is usually taken to be a free parameter controlling the coordinate space scale. NEXT-TO-LEADING ORDER BALITSKY-KOVCHEGOV EQUATION …PHYS. REV. D 102, 074027 (2020) 074027-3 purposes of this paper, the actual form of the initial condition is not relevant. We compare the finite-Ncversion of the BK equation presented above to the equation obtained in the large-Nc limit, which has been studied numerically in Refs. [62,64]. In this limit, one can drop operators suppressed by 1=Nc and the correlators in Eq. (2) become hD1i⟶ Nc→∞hD1iNc→∞¼hSð2Þ x;zihSð2Þ z;yi−hSð2Þ x;yi;ð15Þ hD2;1i⟶ Nc→∞hD2iNc→∞;ð16Þ hD2;2i⟶ Nc→∞hD2iNc→∞;ð17Þ with hD2iNc→∞¼hSð2Þ x;zihSð2Þ z;z0ihSð2Þ z0;yi−hSð2Þ x;zihSð2Þ z;yi:ð18Þ III. SIX-POINT FUNCTIONS IN THE GAUSSIAN APPROXIMATION A. The Gaussian approximation At large Nc, all Wilson line operators present in the NLO BK equation can be expressed solely in terms of dipole correlators, as can be seen from Eqs. (15) and (16). At finite Nc, on the other hand, the higher-point functions hSð2Þ x;zSð2Þ z;yi, hSð2Þ x;zSð2Þ z;z0Sð2Þ z0;yiand hSð6Þ x;z;z0;y;z;z0iare needed. An expression for the four-point function in terms of two-point functions has been derived using the Gaussian approximation (see e.g., [55]). This makes it possible to obtain a closed form for the LO BK equation at finite Nc. The purpose of this work is to compute also the six-point functions using the Gaussian approximation, in order to obtain a closed finiteNcBK equation at NLO accuracy. In the Gaussian approximation, all correlators are parametrized by a single two-point function, and all higherpoint functions can then be expressed in terms of this function. The initial condition for the small-xevolution is usually assumed to be Gaussian (e.g., as in the MV model), but it is not clear a priori that the Gaussian approximation is valid after the evolution. However, numerical studies of the JIMWLK equation [70,71] have not found any indication of major effects breaking the validity of this approximation. We use the diagrammatic notation of Refs. [53,55,72],in which Wilson lines are denoted as ð19Þ ð20Þ The projectile transverse coordinate is x, the lightcone time axis runs from right to left and the blue vertical line represents the target background field. In the Gaussian approximation, the correlator for some Wilson line operator O½Uis approximated as an integral over a parametrization rapidity ηof a single two-point correlator Gu1;u2: hO½Uiη¼exp −1 2Zηd˜ ηZu1;u2 Gu1;u2ð˜ ηÞLa u1La u2O½U: ð21Þ The transverse integrals are denoted as Ru¼Rd2uand La u is a Lie derivative that acts on Wilson lines according to La uUx¼−igδð2Þðx−uÞtaUxð22Þ ð23Þ The structure as an exponential of a two-point function (as in the MV model [67]) is what makes this a “Gaussian” approximation. In practice, for gauge invariant (color singlet) operators, the two-point function Gu1;u2always appears in the linear combination Gu1;u2ðηÞ≔Zηd˜ ηGu1;u2ð˜ ηÞ−1 2ðGu1;u1ð˜ ηÞþGu2;u2ð˜ ηÞÞ: ð24Þ For the integrand, we use the notation G0≔∂ηG. Physical observables only depend on the integrated Gand not on the integrand G0. Thus, for the purpose of our calculation, where we need to relate higher-point functions of Wilson lines to the two-point function, there is some freedom in choosing the parametrization rapidity. We use this freedom in such a way that the ηand transverse coordinate dependences of G0factorize,2as is usually done when employing the Gaussian approximation (see e.g., [54,56,58,70,74–76]) in phenomenological applications. We henceforth omit the explicit ηdependence of Gfor brevity. As an example of the procedure for finding the parametric equation for a correlator using Eq. (21), we illustrate 2This assumption has the effect that the transition matrices MðηÞ[introduced below in Eq. (30)] at different rapidities η commute with each other. This turns the path ordered exponential of MðηÞinto a normal exponential. This, in turn, makes it possible to relate higher-point functions to the two-point function without any further assumptions about the ηdependence of G0.In the terminology of Ref. [55], we use a “rigid exponentiation” instead of the “Gaussian truncation.”The “Gaussian truncation” would imply equating the parametrization rapidity ηwith the evolution rapidity Y. See the related discussion in Refs. [72,73]. T. LAPPI, H. MÄNTYSAARI, and A. RAMNATH PHYS. REV. D 102, 074027 (2020) 074027-4 the steps for the dipole operator hSð2Þ x;yi. A more practical form of Eq. (21) is to write the integral over the parametrization rapidity ηin differential form. Then, we have ∂ηhSð2Þ x;yi¼−MðηÞhSð2Þ x;yi;ð25Þ where Mis the so-called transition matrix formed by operating with the argument of the exponential in Eq. (21) on the operator of interest. In this case, there are four contributions from acting with the Lie derivatives on a product of two Wilson lines: ð26Þ ð27Þ From this sum, we need to factorize out the original operator Ux⊗U† y, so we use the Fierz identity ð28Þ There is only one way to join the endpoints of the Wilson lines into a singlet operator. This is to trace over them, i.e., wedging them between and Doing so and performing the remaining operations in the exponent of Eq. (21), we get ð29Þ where CF¼N2 c−1 2Nc. Since the operators and were normalized, the initial condition at η¼0for the differential equation (25) is given by trivial Wilson lines equal to the identity matrix in the absence of a color field. This is the well-known parametric equation for the dipole correlator [55] in the Gaussian approximation. In the case of n-point correlators larger than the dipole, the operator hO½Ui in Eq. (21) is actually an n×nmatrix of correlators, denoted AðηÞ, and Eq. (25) becomes an n×nmatrix differential equation ∂ηAðηÞ¼−MðηÞAðηÞ:ð30Þ By construction, Mis a symmetric matrix, so there are at most Pn i¼1i¼nðnþ1Þ=2distinct elements, not n×n. For example, a product of six Wilson lines is represented in this notation as ð31Þ (the haphazard assignment of coordinate labels is convenient for the NLO BK equation and will become clear when constructing the transition matrix). The notation here means that this product is actually a matrix with six open indices on the left and another six on the right; we denote them as Since only singlet states are gauge-invariant, the only operators of interest here are singlets. There are six possible ways to join the endpoints of these Wilson lines to form a singlet. Equation (30) is therefore a six-by-six matrix differential equation, as opposed to the much simpler one-dimensional problem illustrated for the dipole correlator. In analogy to the procedure for the dipole operator, the procedure to use Eq. (30) to find parametric equations for the six-point correlators is as follows: (1) Choose a multiplet basis, represented as a column vector B,ofn¼6color structures for the space of all six-point correlators. Each element will have six open color indices, which can contract with the open indices on the left of Uz⊗U† z0⊗Uv⊗ U† y⊗Ux⊗U† w. For example, one choice for an element of Bcould be ð32Þ and the corresponding element for the other end of the Wilson line is then NEXT-TO-LEADING ORDER BALITSKY-KOVCHEGOV EQUATION …PHYS. REV. D 102, 074027 (2020) 074027-5 ð33Þ The prefactor is a normalization constant found by squaring the basis element. (2) Construct the correlator matrix Aby taking BðUz⊗U† z0⊗Uv⊗U† y⊗Ux⊗U† wÞBT. The elements of Bon the left are contracted with the open color indices on the left of the Wilson lines, and the elements of BTwith the open indices on the right. For example, using the basis element shown above, we have for one of the 36 elements in A, ð34Þ (3) Construct the transition matrix Mby summing (for each element in A) all possible one-gluon diagrams obtained with the double Lie derivative operator and rewriting the result in terms of elements of A. (4) Solve Eq. (30) by exponentiating Mto find expressions for each element in A, using as an initial condition the correlator matrix Acorresponding to Wilson lines equal to the identity matrix. B. Choosing a basis Starting from a product of six Wilson lines, there are six ways to form multiplets by joining endpoints in all possible ways: ð35Þ The simplest way to construct an orthonormal basis from these would be to use color algebra to choose ð36Þ The blue lines denote gluons and the last two elements of Brepresent the antisymmetric and symmetric structure constants, respectively, fabc ¼−2itrð½ta;t b;t cÞand dabc ¼2trðfta;t bg;t cÞ. The color factors are dA¼N2 c− 1and Cd¼N2 c−4 Nc. The next step would be to use this basis to construct the correlator matrix Aand the transition matrix M. However, doing so results in a matrix differential equation ∂ηAðηÞ¼−MðηÞAðηÞ, whose complicated solution is the matrix exponential of a six-by-six matrix M. For our case, a better way to proceed is to exploit the structure of the six-point correlators that are actually needed for the NLO BK equation. Since there are only four distinct coordinates in these particular correlators, we make the coordinate assignments ð37Þ It is easy to see from this that there is one way to join the endpoints such that in the limit v→z0,w→z, four Wilson lines cancel (due to unitarity). The result simplifies to a single trace: T. LAPPI, H. MÄNTYSAARI, and A. RAMNATH PHYS. REV. D 102, 074027 (2020) 074027-6 ð38Þ So choosing as one of our basis elements allows one dimension of our six-dimensional space of operators to decouple, giving the equation for the dipole correlator. Similarly, the choice of two more particular basis elements results in two correlators that reduce to four-point functions; one due to the limit v→z0and the other due to the limit w→z. Thus, we can expect to choose a further two basis elements such that two more dimensions decouple from the remaining five, corresponding to the equation for the four-point operators. These two basis elements can be chosen as and We will choose the remaining three basis elements such that they are orthonormal to the three already chosen, resulting in the final basis vector ð39Þ Since this basis is orthonormal, the correlator matrix at the initial condition Aðη¼0Þwill just be the identity matrix. C. Constructing the correlator matrix and the transition matrix Due to this choice of basis ˜ B, the full matrix differential equation (30) now decouples into three independent equations. This allows us to forego exponentiating a six-by-six matrix; at most we will need to exponentiate a three-bythree matrix, which can be done analytically. To form the correlator matrix A, we take the product ˜ BðUz⊗U† z0⊗Uv⊗U† y⊗Ux⊗U† wÞ ˜ BTand set w→z and v→z0. To form the transition matrix, we act with the argument of the exponential in Eq. (21) on the operator Uz⊗U† z0⊗Uv⊗U† y⊗Ux⊗U† w, then wedge the result between the basis vectors and set w→zand v→z0: ˜ B−1 2Zηdη0Zu1;u2 Gu1;u2ðη0ÞLa u1La u2 ×ðUz⊗U† z0⊗Uv⊗U† y⊗Ux⊗U† wÞ ˜ BTw→z v→z0ð40Þ Diagrammatically, this is equivalent to summing all possible ways of attaching one gluon line on the NEXT-TO-LEADING ORDER BALITSKY-KOVCHEGOV EQUATION …PHYS. REV. D 102, 074027 (2020) 074027-7 operator , using the Fierz identity to replace the gluon vertices, and finally closing the Wilson line endpoints on the left and right using the basis vector. For example, the element (6, 6) of the correlator matrix Ais ð41Þ We then sum all the diagrams in which a gluon is attached to this diagram so that it joins any two of the six Wilson lines on the right of the target interaction. Using the Fierz identity in Eq. (28), we may write the resulting expression in terms of diagrams with no gluons. After making the substitutions w→zand v→z0, the result will be a linear combination of elements of the operator matrix AðηÞ. From this linear combination, one can read off the elements of column 6 of the transition matrix M. The explicit expressions for the elements of AðηÞin terms of the Wilson line correlators are shown in Appendix. Performing this procedure for each diagram in AðηÞ,we get the full transition matrix MðηÞjw→z v→z0¼0 B @ M300 0M20 00M1 1 C AðηÞ;ð42Þ where the subscripts refer to the dimension of the submatrix. The first (one-dimensional) transition submatrix is M1ðηÞ¼CFGx;y 0ð43Þ and upon exponentiation, gives the parametric equation for the dipole correlator as shown in Eq. (29). When inverted, this equation can be used to express the two-point function Gx;yin terms of the dipole correlator. This will be needed to evaluate the higher-point functions in terms of the dipole hSð2Þ x;yi. The second (two-dimensional) transition submatrix is M2ðηÞ¼Nc 4M2ð1;1ÞM2ð1;2Þ M2ð1;2ÞM2ð1;1ÞðηÞ;ð44Þ where Mð1;1Þ 2ðηÞ≔G0 x;zþG0 y;z−2 N2 c G0 x;yþG0 x;z0þG0 y;z0;ð45Þ Mð1;2Þ 2ðηÞ≔G0 x;zþG0 y;z−G0 x;z0−G0 y;z0:ð46Þ The matrix differential equation ∂ηA2ðηÞ¼−M2ðηÞA2ðηÞð47Þ then gives a coupled system of 2×2differential equations, out of which 2 are linearly independent, corresponding to the fact that the same transition matrix operates separately on each of the two columns of A2. The exponential solution for this system of equations gives the known parametrization for the four-point correlator with one repeated coordinate [55] hSð2Þ x;zSð2Þ z;yi¼ 1 N2 c e−CFGx;yþ2CF Nc e−CFGx;ye−Nc 2ðGx;zþGy;z−Gx;yÞ: ð48Þ The third and final (three-dimensional) transition submatrix is M3ðηÞ¼0 B B @ Nc 4Γ0 1ffiffiffiffiffiffiffiffi NcCd p 4Γ0 20 ffiffiffiffiffiffiffiffi NcCd p 4Γ0 2 Nc 4Γ0 1−1 ffiffi2 pΓ0 2 0−1 ffiffi2 pΓ0 2Γ0 0 1 C C A ;ð49Þ where Γ0≔CFGx;yþNcGz;z0;ð50Þ Γ1≔Gx;zþGy;z−2 N2 c Gx;yþGx;z0þGy;z0þ2Gz;z0;ð51Þ Γ2≔Gx;z−Gy;z−Gx;z0þGy;z0ð52Þ and the primes on the Γ’s in Eq. (49) denote derivatives in η. Exponentiating this matrix is the last step required to get expressions for the remaining six-point correlators in A3. D. Exponentiating the transition matrix M3 In order to obtain the six-point functions, it is necessary to solve the differential equation ∂ηA3ðηÞ¼−M3ðηÞA3ðηÞ:ð53Þ Solving this equation is equivalent to exponentiating the matrix M3, as shown above in the cases of the twoand four-point functions. To exponentiate M3, we consider two different cases: Γ0 2¼0and Γ0 2≠0. The reason for this will become clear shortly. When Γ0 2¼0,M3in Eq. (49) becomes diagonal and we directly obtain A3ðηÞ¼0 B @ e−Nc 4Γ100 0e−Nc 4Γ10 00e−Γ0 1 C A:ð54Þ When Γ0 2≠0, the matrix elements of A3are calculated by matrix-exponentiating the full M3in Eq. (49), giving T. LAPPI, H. MÄNTYSAARI, and A. RAMNATH PHYS. REV. D 102, 074027 (2020) 074027-8 ðA2Þ The two-dimensional submatrix is ðA3Þ where ðA4Þ ðA5Þ The three-dimensional submatrix is ðA6Þ where ðA7Þ ðA8Þ ðA9Þ NEXT-TO-LEADING ORDER BALITSKY-KOVCHEGOV EQUATION …PHYS. REV. D 102, 074027 (2020) 074027-15 ðA10Þ ðA11Þ ðA12Þ [1] H. Abramowicz et al. (H1 and ZEUS Collaboration), Combination of measurements of inclusive deep inelastic epscattering cross sections and QCD analysis of HERA data, Eur. Phys. J. C 75, 580 (2015). [2] F. Gelis, E. Iancu, J. Jalilian-Marian, and R. Venugopalan, The color glass condensate, Annu. Rev. Nucl. Part. Sci. 60, 463 (2010). [3] J. Jalilian-Marian, A. Kovner, L. D. McLerran, and H. Weigert, The Intrinsic glue distribution at very small x, Phys. Rev. D 55, 5414 (1997). [4] J. Jalilian-Marian, A. Kovner, A. Leonidov, and H. Weigert, The BFKL equation from the Wilson renormalization group, Nucl. Phys. B504, 415 (1997). [5] J. Jalilian-Marian, A. Kovner, A. Leonidov, and H. Weigert, The Wilson renormalization group for low x physics: Towards the high density regime, Phys. Rev. D 59, 014014 (1998). [6] E. Iancu and L. D. McLerran, Saturation and universality in QCD at small x, Phys. Lett. B 510, 145 (2001). [7] Y. V. Kovchegov, Small x F(2) structure function of a nucleus including multiple pomeron exchanges, Phys. Rev. D 60, 034008 (1999). [8] I. Balitsky, Operator expansion for high-energy scattering, Nucl. Phys. B463, 99 (1996). [9] Y. V. Kovchegov and H. Weigert, Triumvirate of running couplings in small-x evolution, Nucl. Phys. A784, 188 (2007). [10] E. Gardi, J. Kuokkanen, K. Rummukainen, and H. Weigert, Running coupling and power corrections in nonlinear evolution at the high-energy limit, Nucl. Phys. A784, 282 (2007). [11] J. L. Albacete and Y. V. Kovchegov, Solving high energy evolution equation including running coupling corrections, Phys. Rev. D 75, 125021 (2007). [12] I. Balitsky, Quark contribution to the small-x evolution of color dipole, Phys. Rev. D 75, 014001 (2007). [13] T. Lappi and H. Mäntysaari, On the running coupling in the JIMWLK equation, Eur. Phys. J. C 73, 2307 (2013). [14] H. Abramowicz et al. (H1 and ZEUS Collaborations), Combination and QCD analysis of charm and beauty production cross-section measurements in deep inelastic ep scattering at HERA, Eur. Phys. J. C 78, 473 (2018). [15] J. L. Albacete, N. Armesto, J. G. Milhano, P. Quiroga-Arias, and C. A. Salgado, AAMQS: A non-linear QCD analysis of new HERA data at small-x including heavy quarks, Eur. Phys. J. C 71, 1705 (2011). [16] J. L. Albacete, J. G. Milhano, P. Quiroga-Arias, and J. Rojo, Linear vs Non-linear QCD evolution: From HERA data to LHC phenomenology, Eur. Phys. J. C 72, 2131 (2012). [17] T. Lappi and H. Mäntysaari, Single inclusive particle production at high energy from HERA data to protonnucleus collisions, Phys. Rev. D 88, 114020 (2013). [18] H. Mäntysaari and B. Schenke, Confronting impact parameter dependent JIMWLK evolution with HERA data, Phys. Rev. D 98, 034013 (2018). [19] A. H. Rezaeian, M. Siddikov, M. Van de Klundert, and R. Venugopalan, Analysis of combined HERA data in the impact-parameter dependent saturation model, Phys. Rev. D 87, 034002 (2013). T. LAPPI, H. MÄNTYSAARI, and A. RAMNATH PHYS. REV. D 102, 074027 (2020) 074027-16 [20] H. Mäntysaari and P. Zurita, In depth analysis of the combined HERA data in the dipole models with and without saturation, Phys. Rev. D 98, 036002 (2018). [21] P. Tribedy and R. Venugopalan, QCD saturation at the LHC: Comparisons of models to p þp and A þA data and predictions for p þPb collisions, Phys. Lett. B 710, 125 (2012); Erratum, Phys. Lett. B 718, 1154 (2013). [22] V. P. Goncalves and M. L. L. da Silva, Probing the Color Glass Condensate in pp collisions at forward rapidities and very low transverse momenta, Nucl. Phys. A906, 28 (2013). [23] B. Duclou´e, T. Lappi, and H. Mäntysaari, Forward J=ψ production in proton-nucleus collisions at high energy, Phys. Rev. D 91, 114005 (2015). [24] B. Duclou´e, T. Lappi, and H. Mäntysaari, Forward J=ψ production at high energy: Centrality dependence and mean transverse momentum, Phys. Rev. D 94, 074031 (2016). [25] J. L. Albacete, P. Guerrero Rodríguez, and Y. Nara, Ultraforward particle production from color glass condensate and Lund fragmentation, Phys. Rev. D 94, 054004 (2016). [26] B. Duclou´e, T. Lappi, and H. Mäntysaari, Isolated photon production in proton-nucleus collisions at forward rapidity, Phys. Rev. D 97, 054023 (2018). [27] H. Mäntysaari and H. Paukkunen, Saturation and forward jets in proton-lead collisions at the LHC, Phys. Rev. D 100, 114029 (2019). [28] E. C. Aschenauer, S. Fazio, J. H. Lee, H. Mantysaari, B. S. Page, B. Schenke, T. Ullrich, R. Venugopalan, and P. Zurita, The electron–ion collider: Assessing the energy dependence of key measurements, Rep. Prog. Phys. 82, 024301 (2019). [29] A. Accardi et al., Electron ion collider: The next QCD frontier - understanding the glue that binds us all, Eur. Phys. J. A 52, 268 (2016). [30] J. Abelleira Fernandez et al. (LHeC Study Group Collaboration), A large hadron electron collider at CERN: Report on the physics and design concepts for machine and detector, J. Phys. G 39, 075001 (2012). [31] X. Chen, A plan for electron ion collider in China, Proc. Sci. DIS2018 (2018) 170 [arXiv:1809.00448]. [32] C. A. Bertulani, S. R. Klein, and J. Nystrand, Physics of ultra-peripheral nuclear collisions, Annu. Rev. Nucl. Part. Sci. 55, 271 (2005). [33] S. R. Klein and H. Mäntysaari, Imaging the nucleus with high-energy photons, Nat. Rev. Phys. 1, 662 (2019). [34] I. Balitsky and G. A. Chirilli, Next-to-leading order evolution of color dipoles, Phys. Rev. D 77, 014019 (2008). [35] I. Balitsky and G. A. Chirilli, Rapidity evolution of Wilson lines at the next-to-leading order, Phys. Rev. D 88, 111501 (2013). [36] A. Kovner, M. Lublinsky, and Y. Mulian, Jalilian-Marian, Iancu, McLerran, Weigert, Leonidov, Kovner evolution at next to leading order, Phys. Rev. D 89, 061704 (2014). [37] B. Duclou´e, H. Hänninen, T. Lappi, and Y. Zhu, Deep inelastic scattering in the dipole picture at next-to-leading order, Phys. Rev. D 96, 094017 (2017). [38] G. Beuf, Dipole factorization for DIS at NLO: Combining the q¯ qand q¯ qg contributions, Phys. Rev. D 96, 074033 (2017). [39] G. Beuf, Dipole factorization for DIS at NLO: Loop correction to the γ T;L →q¯ qlight-front wave functions, Phys. Rev. D 94, 054016 (2016). [40] I. Balitsky and G. A. Chirilli, Photon impact factor in the next-to-leading order, Phys. Rev. D 83, 031502 (2011). [41] H. Hänninen, T. Lappi, and R. Paatelainen, One-loop corrections to light cone wave functions: The dipole picture DIS cross section, Ann. Phys. (Amsterdam) 393, 358 (2018). [42] R. Boussarie, A. V. Grabovsky, D. Yu. Ivanov, L. Szymanowski, and S. Wallon, Next-to-Leading Order Computation of Exclusive Diffractive Light Vector Meson Production in a Saturation Framework, Phys. Rev. Lett. 119, 072002 (2017). [43] M. Escobedo and T. Lappi, Dipole picture and the nonrelativistic expansion, Phys. Rev. D 101, 034030 (2020). [44] T. Lappi, H. Mäntysaari, and J. Penttala, Relativistic corrections to the vector meson light front wave function, arXiv:2006.02830. [45] G. A. Chirilli, B.-W. Xiao, and F. Yuan, Inclusive hadron productions in pA collisions, Phys. Rev. D 86, 054005 (2012). [46] B. Duclou´e, E. Iancu, T. Lappi, A. H. Mueller, G. Soyez, D. N. Triantafyllopoulos, and Y. Zhu, Use of a running coupling in the NLO calculation of forward hadron production, Phys. Rev. D 97, 054020 (2018). [47] B. Duclou´e, T. Lappi, and Y. Zhu, Single inclusive forward hadron production at next-to-leading order, Phys. Rev. D 93, 114016 (2016). [48] K. Watanabe, B.-W. Xiao, F. Yuan, and D. Zaslavsky, Implementing the exact kinematical constraint in the saturation formalism, Phys. Rev. D 92, 034026 (2015). [49] T. Altinoluk, N. Armesto, G. Beuf, A. Kovner, and M. Lublinsky, Single-inclusive particle production in protonnucleus collisions at next-to-leading order in the hybrid formalism, Phys. Rev. D 91, 094016 (2015). [50] A. M. Stasto, B.-W. Xiao, and D. Zaslavsky, Towards the Test of Saturation Physics Beyond Leading Logarithm, Phys. Rev. Lett. 112, 012302 (2014). [51] H.-Y. Liu, Y.-Q. Ma, and K.-T. Chao, Improvement for color glass condensate factorization: Single hadron production in pA collisions at next-to-leading order, Phys. Rev. D 100, 071503 (2019). [52] K. Rummukainen and H. Weigert, Universal features of JIMWLK and BK evolution at small x, Nucl. Phys. A739, 183 (2004). [53] Y. V. Kovchegov, J. Kuokkanen, K. Rummukainen, and H. Weigert, Subleading-N(c) corrections in non-linear small-x evolution, Nucl. Phys. A823, 47 (2009). [54] H. Fujii, F. Gelis, and R. Venugopalan, Quark pair production in high energy pA collisions: General features, Nucl. Phys. A780, 146 (2006). [55] C. Marquet and H. Weigert, New observables to test the color glass condensate beyond the large-Nclimit, Nucl. Phys. A843, 68 (2010). [56] F. Dominguez, C. Marquet, B.-W. Xiao, and F. Yuan, Universality of unintegrated gluon distributions at small x, Phys. Rev. D 83, 105005 (2011). NEXT-TO-LEADING ORDER BALITSKY-KOVCHEGOV EQUATION …PHYS. REV. D 102, 074027 (2020) 074027-17 [57] C. Marquet, Forward inclusive dijet production and azimuthal correlations in p(A) collisions, Nucl. Phys. A796,41 (2007). [58] T. Lappi and H. Mantysaari, Forward dihadron correlations in deuteron-gold collisions with the Gaussian approximation of JIMWLK, Nucl. Phys. A908, 51 (2013). [59] K. Dusling, M. Mace, and R. Venugopalan, Parton model description of multiparticle azimuthal correlations in pA collisions, Phys. Rev. D 97, 016014 (2018). [60] F. Dominguez, C. Marquet, A. M. Stasto, and B.-W. Xiao, Universality of multiparticle production in QCD at high energies, Phys. Rev. D 87, 034007 (2013). [61] B. Duclou´e, E. Iancu, A. H. Mueller, G. Soyez, and D. N. Triantafyllopoulos, Non-linear evolution in QCD at highenergy beyond leading order, J. High Energy Phys. 04 (2019) 081. [62] T. Lappi and H. Mäntysaari, Direct numerical solution of the coordinate space Balitsky-Kovchegov equation at next to leading order, Phys. Rev. D 91, 074016 (2015). [63] E. Iancu, J. D. Madrigal, A. H. Mueller, G. Soyez, and D. N. Triantafyllopoulos, Resumming double logarithms in the QCD evolution of color dipoles, Phys. Lett. B 744, 293 (2015). [64] T. Lappi and H. Mäntysaari, Next-to-leading order BalitskyKovchegov equation with resummation, Phys. Rev. D 93, 094004 (2016). [65] G. Beuf, Improving the kinematics for low-xQCD evolution equations in coordinate space, Phys. Rev. D 89, 074039 (2014). [66] E. Iancu, J. D. Madrigal, A. H. Mueller, G. Soyez, and D. N. Triantafyllopoulos, Collinearly-improved BK evolution meets the HERA data, Phys. Lett. B 750, 643 (2015). [67] L. D. McLerran and R. Venugopalan, Computing quark and gluon distribution functions for very large nuclei, Phys. Rev. D49, 2233 (1994). [68] L. D. McLerran and R. Venugopalan, Boost covariant gluon distributions in large nuclei, Phys. Lett. B 424,15 (1998). [69] K. J. Golec-Biernat and M. Wusthoff, Saturation effects in deep inelastic scattering at low Q2and its implications on diffraction, Phys. Rev. D 59, 014017 (1998). [70] A. Dumitru, J. Jalilian-Marian, T. Lappi, B. Schenke, and R. Venugopalan, Renormalization group evolution of multigluon correlators in high energy QCD, Phys. Lett. B 706, 219 (2011). [71] T. Lappi, B. Schenke, S. Schlichting, and R. Venugopalan, Tracing the origin of azimuthal gluon correlations in the color glass condensate, J. High Energy Phys. 01 (2016) 061. [72] T. Lappi, A. Ramnath, K. Rummukainen, and H. Weigert, JIMWLK evolution of the odderon, Phys. Rev. D 94, 054014 (2016). [73] E. Iancu and D. Triantafyllopoulos, Higher-point correlations from the JIMWLK evolution, J. High Energy Phys. 11 (2011) 105. [74] J. P. Blaizot, F. Gelis, and R. Venugopalan, High-energy pA collisions in the color glass condensate approach. 1. Gluon production and the Cronin effect, Nucl. Phys. A743,13 (2004). [75] J. P. Blaizot, F. Gelis, and R. Venugopalan, High-energy pA collisions in the color glass condensate approach. 2. Quark production, Nucl. Phys. A743, 57 (2004). [76] F. Dominguez, C. Marquet, and B. Wu, On multiple scatterings of mesons in hot and cold QCD matter, Nucl. Phys. A823, 99 (2009). T. LAPPI, H. MÄNTYSAARI, and A. RAMNATH PHYS. REV. D 102, 074027 (2020) 074027-18