$\gamma^* \gamma^*$ Cross Section at NLO and Properties of the BFKL Evolution at Higher Orders

We obtain a simple analytic expression for the high energy $\gamma^* \gamma^*$ scattering cross section at the next-to-leading order in the logarithms-of-energy power counting. To this end we employ the eigenfunctions of the NLO BFKL equation constructed in our previous paper. We also construct the eigenfunctions of the NNLO BFKL kernel and obtain a general form of the solution for the NNLO BFKL equation, which confirms the ansatz proposed in our previous paper.


I. INTRODUCTION
Understanding the high energy dynamics of strong interactions is important both for the description of the strong interaction data reported by the current and future accelerators worldwide, and to improve our theoretical understanding of quantum chromodynamics (QCD). Significant progress towards achieving these goals has been accomplished in the recent decades by the physics of saturation and Color Glass Condensate (CGC) (see [1][2][3][4][5][6] for reviews). This approach is based on the existence of an intrinsic hard momentum scale which characterizes hadronic and nuclear wave functions at high energy -the saturation scale Q s [7]. At small values of Bjorken x and/or for large nuclei the saturation scale Q s is much larger than the QCD confinement scale, justifying the use of the small-coupling expansion since α s (Q 2 s ) ≪ 1, and allowing for first-principles calculations of many high energy scattering cross sections and of other related observables.
Calculations of observables in the saturation/CGC framework consist of two steps: (i) first one has to calculate the observables in the quasi-classical Glauber-Mueller (GM) [8] / McLerran-Venugopalan (MV) approximation [9][10][11] and then (ii) evolve the result with energy s using the non-linear Balitsky-Kovchegov (BK) [12][13][14][15] or Jalilian-Marian-Iancu-McLerran-Weigert-Leonidov-Kovner (JIMWLK) [16][17][18][19] evolution equations. It appears that to successfully describe and, sometimes, predict a host of deep inelastic scattering (DIS), hadronic and nuclear collisions data (see e.g. [20][21][22]) one needs to include the running-coupling corrections [23][24][25][26][27] into the leading-ln s BK and JIMWLK evolution equations. (Henceforth we will refer to the leading-ln s BK and JIMWLK equations either as leading-order (LO) BK and JIMWLK equations.) The next-to-leading order (NLO) corrections to both the BK [28] and JIMWLK [29,30] evolution equations are known, but have not yet been implemented in phenomenological applications. At the same time keeping the higherorder corrections under theoretical control is essential for improving the precision of the saturation/CGC physics predictions. Perhaps more importantly, an agreement of any theory with the data is established only once the uncertainty due to higher-order corrections is understood and quantified, demonstrating explicitly that higher-order corrections do not significantly affect the existing agreement of the theory with the data. It is, therefore, very important to understand both the structure and magnitude of the higher-order corrections to the small-x evolution equations.
Since no exact analytic solutions for the nonlinear LO BK and JIMWLK equations are known, it is easier to start organizing the higher-order corrections for the linear Balitsky-Fadin-Kuraev-Lipatov (BFKL) [31][32][33] evolution equation, which predated both the BK and JIMWLK equations by over two decades: in the linear (low parton density) regime both the BK and JIMWLK equations reduce to BFKL. The LO BFKL equation's solution is well-known and was constructed almost simultaneously with the equation itself [33]. The kernel of the BFKL equation is known up to the next-to-leading order [34,35]. In the fifteen years since the construction of the NLO BFKL equation numerous efforts have been made to understand the features of its solution: non-Regge behavior was found in the amplitude resulting from the NLO BFKL equation [36][37][38][39]; it was also shown that the NLO correction to the BFKL intercept is large and negative [34,[40][41][42], though higher-order collinear divergences are likely to reduce the size of this correction [43][44][45].
More recently the eigenfunctions of the NLO BFKL kernel were constructed in our previous paper [46] by using a perturbative expansion around the eigenfunctions of the LO BFKL kernel. This method yields a systematic way of constructing any-order BFKL kernel eigenfunctions, therefore allowing one to find the solution of any-order BFKL equation. (See [47] for an alternative, though perhaps related way of solving the NLO BFKL equation.) In the general form of the BFKL solution found in [46] the higher-order corrections enter through the perturbative expansions in the intercept and in the eigenfunctions. This is an essential difference between the BFKL solution in QCD and in a conformal field theory, such as the N = 4 super-Yang-Mills (SYM) theory: in the latter, conformal symmetry uniquely fixes the eigenfunctions of the BFKL kernel to be the eigenfunctions of the Casimir operators of the Möbius group, E n,ν [48], such that the higher-order corrections enter the SYM BFKL solution only through perturbative expansion in the intercept [49][50][51][52]. Note also a similar difference from the solution of the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) evolution equation [53][54][55], where the higher-order corrections in the strong coupling α s only enter in the anomalous dimensions of the quark and gluon distribution operators [55][56][57][58][59][60]. The eigenfunctions found in [46] allowed us to construct a simple form of the NLO BFKL Green function and to conjecture an ansatz for the next-to-next-to-leading order (NNLO) BFKL Green function in the same reference.
The aim of the present paper is twofold. One goal is to construct the NNLO BFKL kernel eigenfunctions in QCD by finding perturbative corrections to the NLO BFKL kernel eigenfunctions of [46]. This would allow us to verify the NNLO BFKL Green function ansatz proposed in [46]. The second goal is to build on the existing work [52,61,62] along with, of course, [46], to construct the cross section for high-energy γ * γ * scattering at NLO. We use the standard BFKL power counting in which where Y ∼ ln s is the rapidity variable and α s is the strong coupling constant. In this power counting the LO γ * γ * scattering mediated by the LO BFKL exchange is O(α 2 s ), while the NLO correction we construct below is O(α 3 s ). The resulting NLO γ * γ * scattering cross section would not only allow for a comparison of the prediction with the existing and future experimental data from linear colliders, but would also allow one to see the structure of higher-order corrections to the γ * γ * scattering by obtaining a general form of this cross section (in the linear regime).
The outline of the paper is as follows: we begin in Sec. II by generalizing the NLO BFKL kernel Green functions to the azimuthal-angle dependent case, not considered originally in [46]. This generalization is rather straightforward for a reader familiar with the technique of [46], and results in the eigenfunctions H n,ν (k) given in Eq. (13). We continue in Sec. III where we construct the NLO γ * γ * scattering cross section. The answer is particularly simple and is given in Eqs. (58) and (59) as a convolution of the LO+NLO evolution with the LO+NLO impact factors and LO+NLO daughter dipole-dipole scattering. It appears that higher-order corrections, in the linear regime considered, only enter the expression for the γ * γ * cross section either through the intercept of the evolution, the impact factor or the daughter dipole-dipole scattering amplitude. The γ * γ * cross section is plotted in Figs. 3 and 4. Finally, the NNLO BFKL kernel eigenfunctions are constructed via a direct calculation in Sec. IV and are presented in Eq. (93) below (along with Eqs. (68) and (100)). The corresponding eigenvalue is given by Eq. (99). Using the eigenfunctions to construct the NNLO BFKL Green function (105) we reproduce (and, therefore, confirm) the ansatz proposed in [46]. (Indeed our NNLO BFKL Green function (105) depends on the NNLO BFKL intercept χ 2 (ν), which is not known at present.) We present conclusions and outlook in Sec. V.

II. GENERALIZING NLO BFKL EIGENFUNCTIONS TO THE AZIMUTHALLY-DEPENDENT CASE
The NLO BFKL eigenfunction were constructed in [46] in the azimuthal-angle independent approximation. However, those results are easy to generalize to the case with non-trivial azimuthal angular dependence.
Consider an arbitrary-order BFKL equation for the Green function G q ⊥ , k ′ ⊥ , Y . The initial condition is Here Y is the rapidity variable and k ⊥ , k ′ ⊥ and q ⊥ denote transverse momenta. The BFKL kernel can be expanded in the powers of the strong coupling α µ , such that where µ is an arbitrary renormalization scale and N c is the number of colors. To find the NLO BFKL eigenfunctions we need to know the action of the LO+NLO kernel defined by on the LO BFKL eigenfunctions q −1+2 i ν e i n φq , where q = | q ⊥ |, φ q is the azimuthal angle of q ⊥ with respect to some chosen direction in the transverse plane, ν is a real parameter and n is an integer. This projection was calculated in [49,61] yielding where k = | k ⊥ | and φ k is the azimuthal angle of vector k ⊥ with respect to the same chosen direction in the transverse plane. The LO BFKL kernel eigenvalue is with integer n and real ν, and where ψ(z) = d ln Γ(z)/dz as usual. The prime denotes derivatives with respect to ν, such that χ ′ 0 (n, ν) = ∂χ 0 (n, ν)/∂ν. The one-loop running of the strong coupling is defined as where N f is the number of quark flavors. The real part of the projection of the NLO BFKL kernel on the LO BFKL eigenfunctions is given by [49,61] where χ ′′ 0 (n, ν) = ∂ 2 χ 0 (n, ν)/∂ν 2 and Φ(n, ν) To construct the eigenfunctions of K LO+NLO one could follow the steps employed in [46] for the azimuthallysymmetric case. One may notice, however, that for n = 0 the construction of eigenfunctions is identical to the n = 0 case of [46] with the simple replacements Performing the substitutions (12) in Eq. (47) of Ref. [46] we obtain the eigenfunctions of the LO+NLO BFKL kernel in the azimuthal angle-dependent case Similar to [46] one can show that the functions H n,ν ( k ⊥ ) satisfy the completeness and orthogonality relations up to (and including) O(α µ ), which is the accuracy of the NLO expansion for the eigenfunctions. (The asterisk denotes complex conjugation.) The completeness relation (14) is already satisfied by the LO eigenfunctions and is not affected by the perturbative corrections used to construct the NLO BFKL eigenfunctions (13). Therefore, the space of functions upon which the eigenfunctions H n,ν ( k ⊥ ) form a complete basis is the same as the one of the power-like (LO BFKL) eigenfunctions. By analogy with [46] the eigenvalues of K LO+NLO corresponding to the eigenfunctions H n,ν ( k ⊥ ) are such that the Green function in Eq. (2) is This accomplishes the generalization of the NLO BFKL eigenfunctions to the case of non-trivial azimuthal dependence. Below, when discussing the azimuthally-symmetric n = 0 case we will use the simplified notation which also connects to the notation used in [46].
The three main ingredients for the determination of the NLO cross section for γ * γ * scattering are the impact factor (light-cone wave function squared), calculated up to NLO [61], the solution of the NLO BFKL evolution equation [46], and the energy-independent forward scattering amplitude at NLO of the two fundamental "daughter" color-dipoles (quark-antiquark pairs) produced by the two virtual photons through the evolution before the interaction [52,62]. The convolution of these three contributions yields the total inclusive γ * γ * scattering cross section. Fig. 1 provides a diagrammatic representation of the factorization of the γ * γ * scattering cross section.
The Operator Product Expansion (OPE) at high-energy [12,51] provides a natural framework were the scattering amplitude factorizes in these three main contributions. In the next subsection we review the logic of the OPE at highenergy before performing the actual calculation of the LO+NLO γ * γ * scattering cross section in the two subsequent subsections.

A. OPE at high-energy
In this section we review the high-energy OPE [12] which we will employ to calculate the NLO γ * γ * cross section below.
Let us consider the forward scattering amplitude of two virtual photons with four-momenta q 1 and q 2 and polarizations λ 1 and λ 2 represented by four electromagnetic currents The left panel contains a diagrammatic representation of the γ * γ * scattering amplitude factorized into the impact factors, small-x evolution, and the (energy-independent) daughter dipole-dipole scattering. This structure of the factorization is suggested by the high-energy OPE applied to the electromagnetic currents. In the right panel, the structure of the LO+NLO impact factor is illustrated by some of the diagrams in the shock-wave formalism originally calculated at NLO in [61,63]. Examples of diagrams contributing to the daughter dipole-dipole scattering at LO and NLO are also shown in the right panel.
Here ε λ ρ (q) are the gluon polarization vectors, x ± = (x 0 + x 3 )/ √ 2, and j σ are the electromagnetic currents. Since we are interested in the γ * γ * cross section at high energy (in the Regge limit), we factorize (19) in the following way where S(A) is the part of the QCD action which depends only on the gluon field A µ , det(i∇) is the fermionic determinant with ∇ the covariant derivative, angle brackets . . . A denote the expectation value in the background gluon field A µ , and the normalization factor N is In the following we omit the factor f e 2 f 2 due to the electromagnetic charges of the quarks e f of flavor f , which will need to be reinstated for future phenomenological applications of our results.
The factorization of Eq. (20) is simply due to the fact that in the high-energy asymptotics the dominant contribution to the γ * γ * forward scattering amplitude should contain at least a single two-gluon exchange in the t-channel. A typical dominant high-energy diagram that can be obtained from Eq. (20) is shown in Fig. 2a. An energy-suppressed diagram of the type we neglected in writing Eq. (20) is shown in Fig. 2b, and belongs to the class of the so-called "box" diagrams, where there is a quark t-channel exchange spanning the whole rapidity interval of the scattering process. Such diagrams are suppressed by a power of center-of-mass energy squared s compared to the leading graphs in Fig. 2a, and can be neglected in the high-energy asymptotics.
Each term in the angle brackets . . . A in Eq. (20) is evaluated using the OPE at high energy [12] and the background field technique. The coefficient functions of the OPE of two electromagnetic currents at high energy is the photon impact factor (PIF) as originally defined in [33]. The PIF is convoluted with a matrix element of an operator made out of Wilson lines along the light cone.
Since we are interested in the NLO γ * γ * cross section, we need the LO and NLO PIF. The calculation of NLO PIF has been performed in [61,63]. We will only employ the results of that work in our calculation and refer the reader to the original papers [61,63] for the details of the PIF calculation. We stress, however, that at NLO the  (20) which contributes to the high-energy asymptotics of γ * γ * scattering. The right panel shows an example of the "box" diagram that can not be generated from Eq. (20): box diagrams are suppressed at high energies by a power of energy and, therefore, can be neglected. However, as we will see in the numerical results section, such types of diagrams become relevant, and therefore not anymore negligible at low rapidity where energy is not high enough to suppress such contributions. PIF can be contaminated by the rapidity-dependent (or energy-dependent) terms. As explained in [51], when the high-energy OPE is performed in terms of composite Wilson-line operators, by redefining the operator at hand all the energy dependent terms can be shifted into the matrix elements of Wilson lines, thus leaving the NLO impact factor energy-scale invariant. This is what has been done in the PIF calculation of [61,63]. In addition, the composite Wilson lines operator obtained by the shifting procedure of [61,63] renders the impact factor conformally invariant also at NLO: the energy-dependent terms clearly would have broken conformal invariance. The high-energy OPE in terms of composite Wilson line operator provides an operatorial and systematic procedure to factorize scattering amplitudes into coefficient functions which are energy-scale invariant and matrix elements of composite Wilson line operators that encode the energy dependence of the amplitude.
In γ * γ * scattering at high energy, each of the two virtual photons can be treated at the same time as the projectile and as the target due to the symmetry of the process. We will work in a frame in which the photons collide head-on, such that their four-momenta are with q + 1 and q − 2 very large. Further, we will work in a gauge which reduces to the A + = 0 light-cone gauge near the x − = 0 light cone, and to the A − = 0 light-cone gauge near the x + = 0 light cone, that is, it reduces to the light-cone gauges for each of the photons: examples of such gauges include the Coulomb ∇· A = 0 gauge and the A 0 = 0 temporal gauge. In such gauge, each virtual photon, long before the interaction, splits into a quark-antiquark pair which then traverses the gluonic field produced by the target (the other qq pair): this interaction is described by a color dipole made of two fundamental Wilson lines scattering on another fundamental dipole. In DIS the small-x evolution of Wilson lines is described by the Balitsky-JIMWLK non-linear equation [12,[16][17][18][19]; in the case of a fundamental dipole in the large-N c limit the evolution equation corresponds to the Balitsky-Kovchegov equation [12][13][14][15] and its linearization corresponds to the BFKL equation [31][32][33]. While it is not presently clear which evolution equation describes Wilson-line correlators in γ * γ * scattering at high energy, it is clear that in the linear regime (see Sec. IV for the definition of such regime) the evolution of a dipole operator is given by the BFKL equation.
We are interested in the γ * γ * scattering cross section in the linear case where the resummation in the leading logarithmic approximation in the Regge limit is achieved by the linearization of the evolution equation for the color Wilson line operator defined as where U η ( x ⊥ ) is the Wilson line operator in the fundamental representation, which, for a projectile moving along the x + light cone is defined by and x and y are the points of interaction of the quark and anti-quark pair with the gluonic external field of the target; in the center-of-mass frame the gluonic external field of the target reduces to a shock wave due to Lorentz contraction and time dilation. The Wilson line operator depends on the rapidity η, which may be included in Eq. (24) either by changing the slope of the Wilson line to also include the x − direction or by a rigid cut-off imposed on the longitudinal momenta of the gluons when calculating matrix elements of the Wilson line operators. In the DIS case the NLO evolution of the color dipole (23) with respect to rapidity η has been calculated in [64] where its linearization was proven to coincide with the known result for the NLO BFKL equation [34,35]. At NLO the high-energy OPE is performed in terms of the composite operator [U a ( x ⊥ , y ⊥ )] comp which differs from the one in Eq. (23) by the addition of a counter term that restores conformal symmetry present in the leading order impact factor [51,61,63]. Note also, that the evolution parameter for the composite operator is not the rapidity parameter η from Eq. (23), but is the new parameter a which ensures that the evolution equation of the composite operator [U a ( x ⊥ , y ⊥ )] comp with respect to a is the same as the one for U η ( x ⊥ , y ⊥ ) with respect to η at LO. At NLO the evolution equation for [U a ( x ⊥ , y ⊥ )] comp with respect to a has a conformal kernel in N = 4 SYM; in QCD the NLO evolution kernel for [U a ( x ⊥ , y ⊥ )] comp consists of a conformally-invariant piece and a running-coupling term.
For DIS the momentum-space high-energy OPE of two electromagnetic currents (in the linearized approximation, that is, in absence of nonlinear saturation effects) is [61,63] where I ρσ (q, k) is the photon impact factor and the Wilson-line operator in the matrix element is (see Appendix A for details) Note that in Eq. (25) the composite operator [U a ( k ⊥ , z ⊥ )] comp has to be used instead of U η ( k ⊥ , z ⊥ ) at NLO and beyond. The background field in the averaging . . . A is due to the target, which could be a proton or a nucleus. The impact factor I ρσ (q, k) is known up to NLO [61,63]. At the LO and NLO the impact factor can be written as a Mellin transform with [61,63] where the functions F i (ν) with i = 1, 2, 3 are defined in [61], k = | k ⊥ |, and 1 1 Note that P ρσ 3 depends only on the direction of k ⊥ , and not on its magnitude.
with p µ 2 = (0, p − 2 , 0 ⊥ ). Note that the tensor structure P ρσ 3 in (28) will not contribute in the case of total cross section of γ * γ * scattering we are about to investigate. This is due to the fact that scattering of two unpolarized photons in the kinematics of (22) has no preferred transverse direction: averaging P ρσ 3 over the angles of k ρ ⊥ gives zero. This conclusion also holds for the scattering of longitudinally polarized photons and for the scattering of transversely polarized photons summed over transverse polarizations.
We now turn our attention back to γ * γ * scattering. Using Eqs. (25) and (27) in the scattering amplitude (20) yields and a 1 , a 2 are fixed so that the NLO impact factors in coordinate space are conformally invariant [61]; in momentum space a 1 , a 2 are related to the rapidities of the two colliding virtual photons, a 1 = s/Q 2 1 and a 2 = s/Q 2 2 with s = (q 1 + q 2 ) 2 the center-of-mass energy squared. Note that (cf. Eq. (26) along with Eq. (20)) Opening the square brackets in Eq. (30) we can write the γ * γ * forward scattering amplitude as Using the OPE in terms of composite Wilson line operators, the forward scattering amplitude (33) is now factorized into the energy-scale invariant impact factorsĨ ρσ LO+NLO for each of the virtual photons, the evolution of the composite Wilson line operators [U a ] comp convoluted with each impact factor and the scattering amplitude of the composite (31) and (32) which has to be calculated at NLO (without evolution). Notice that, the NLO impact factor calculated in [65][66][67], which is provided as a combination of numerical and analytic results, would differ from the one calculated in [61,63]. This is because the NLO dipole-dipole scattering, which, in the OPE language, is a universal factor independent of the process under consideration, has to be calculated only once, while in [65][66][67] it is included in the definition of the impact factor. Thus, in [65][66][67] the scattering cross-sections are factorized in a different way and the contribution of the dipole-dipole scattering has to be re-calculated as part of the impact factor for every different process considered.
Before proceeding with the NLO calculation, it is instructive to rederive the LO γ * γ * cross section. Such calculation would also allow one to better understand the issues that have prevented the calculation of γ * γ * cross section at NLO until now. The LO γ * γ * cross section is defined as the O(α 2 s ) contribution in the power counting of Eq. (1). Since at LO we do not need to use the composite operator Eq. (33) reduces to with the LO impact factor [61] I ρσ We now need to calculate the dipole-dipole scattering amplitude U η1 ( k 1⊥ , z ⊥ ) U η2 ( k 2⊥ , 0 ⊥ ) at LO including the high-energy leading logarithmic (LO BFKL) resummation. Each operator k 2 i U ηi needs to be evolved using the LO BFKL evolution. To do so we employ the fact that transverse momentum powers are eigenfunctions of the LO BFKL. Using Mellin representation (see Eqs. (A9) and (A11) in Appendix A) with we write the LO BFKL evolution for the (Mellin moment of the) dipole amplitude as where η 0 is some initial rapidity value and χ 0 (ν) is defined in Eqs. (12) and (8).
Using Eqs. (36) and (38) in (34) we have where η 1 − η 2 = ln s Q1Q2 with s = (q 1 + q 2 ) 2 the center-of-mass energy squared. Since each U η (ν) in Eq. (39) contains an integral over impact parameters (see (37)), while in the γ * γ * scattering at hand there is only one integration over the impact parameter (see Eq. (34)), we have to divide out the infinite factor of S ⊥ = d 2 z from the correlator in Eq. (39).
We now need the scattering of two color-dipoles at LO in Mellin space. This calculation is provided in [52] but in N = 4 SYM theory. To obtain the QCD expression we have to modify only the color factors getting Using (40) in (39) and making use of (A9) and (A10) yields where we have used the fact thatĨ ρσ LO (q, −ν) =Ĩ ρσ LO (q, ν) and χ 0 (ν) = χ 0 (−ν). To obtain the total LO γ * γ * cross section we use the optical theorem to get where we have replaced η 1 − η 2 by ln s Q1Q2 . To calculate the NLO γ * γ * cross section one could repeat at the next order inᾱ µ the steps we have just performed to get the LO cross section. At NLO, however, it is well known that the impact factor may get contributions of energy-dependent terms; hence, if we want to preserve the factorization structure of the LO scattering amplitude, we need to define the impact factor in such a way that the NLO impact factor is not dependent on the energy scale. This result is provided by the high-energy OPE in terms of the composite Wilson line operator. The composite Wilson line operator is the dipole operator defined in (23) with counterterms that restore the conformal invariance violated by the energy-dependent terms present at NLO [61]. The NLO impact factor is defined as the coefficient function in front of the composite Wilson line operator. At the moment, the high-energy OPE in terms of composite Wilson line operator [51] represents the only known way to systematically define and calculate the impact factor at any order as the coefficient function in the OPE language. The energy-independent NLO impact factor was found in [61] allowing us to proceed with the calculation of the NLO γ * γ * cross section.
There is another issue which had to be resolved before the calculation of the NLO γ * γ * cross section is possible. One had to solve the NLO BFKL equation in order to also include NLO evolution into the cross section. The NLO BFKL equation has been solved recently in [46] by constructing the NLO BFKL eigenfunctions and the corresponding eigenvalues.
We now have at our hand all the necessary ingredients to proceed to the calculation of the NLO γ * γ * cross section.
C. Assembling the NLO γ * γ * cross section We start with Eq. (33), which, with the help of Eq. (27) we rewrite as where As shown in [61], linearization of the NLO evolution equation for the operator coincides with the NLO BFKL equation obtained in [68]. Therefore, the eigenfunctions (13) are the eigenfunctions of the kernel of the linearized equation for L a (k). Considering only the n = 0 contribution which dominates in high energy scattering, we use the completeness relation (14) for n = 0, to rewrite Eq. (44) as with a projection of the impact factor on the LO+NLO BFKL eigenfunctions.
Since, as we have just mentioned, the functions H ν (k) are the eigenfunctions of the evolution equation for L a (k) [61], the linearized evolution of L a (k) for n = 0 can be written as where and We thus rewrite Eq. (48) as (cf. Eq. (39)) Here (1/2) ln a 0 is an arbitrary rapidity for the LO+NLO dipole-dipole scattering while (1/2) ln(a 1 a 2 ) = ln s Q1Q2 = η 1 − η 2 . Note that the expression (53) for the forward amplitude is quite general, and is likely to be valid at any order in α s within the linear evolution approximation.
We see that the eigenfunctions of LO+NLO BFKL allowed us to write down the evolution effects explicitly in Eq. (53). To complete the calculation, we now need to find the projections I ρσ (q, ν) of the impact factors onto the eigenfunctions H ν (k) (see Eq. (49)) and to construct the correlator L η0 (ν 1 ) L η0 (ν 2 ) at the LO+NLO level. These quantities are calculated in the Appendix B. The result for the dipole-dipole forward scattering amplitude is where the NLO correction to the eikonal dipole-dipole scattering was calculated in [52] for N = 4 SYM, which we modified to work for QCD obtaining 2 Let us stress here, that it is important to include the NLO correction to the "daughter" dipole-dipole scattering to obtain a complete expression for the γ * γ * cross section. Note the difference between the NLO dipole-dipole scattering and the NLO corrections to the impact factor: the latter include the O(α s ) correction to the γ * → qq light-cone wave function squared due to non-eikonal q → qG splittings and mergers; the former includes O(α 3 s ) corrections to the leading-order (O(α 2 s )) scattering of two dipoles made out of pairs of eikonal Wilson lines. (Note that, while we put rapidities of the two operators in Eq. (54) to be equal, this is done in order to exclude small-x evolution corrections; the "daughter" dipole-dipole scattering is still assumed to be high-energy.) The calculation in Appendix B for the LO+NLO impact factor projection yields whereĨ ρσ LO+NLO (q, ν) is the Mellin transform of the LO+NLO impact factor defined in Eq. (27) and ∂ ν ≡ ∂/∂ν. Substituting Eqs. (54) and (56) into Eq. (53) we obtain (see Appendix B for details) where χ 0 (ν) = χ 0 (0, ν) and χ 1 (ν) = χ 1 (0, ν). We have also made use of the one-loop running coupling α s ( with the NLO precision. Note that the term in the last line of Eq. (57) does not contribute if the polarizations of the two virtual photons are either identical or summed over.
The prefactor of Eq. (57) contains α s (Q 2 1 ) α s (Q 2 2 ): with the NLO precision this term is indistinguishable from α 2 s (Q 1 Q 2 ). However, we made this choice because it is known that in the Q 1 ≫ Q 2 regime only one of the couplings should become very small while the other one would encode the non-perturbative features of the "target" like in DIS. An explicit higher-order calculation (see Appendix A of [70]) employing the Brodsky-Lepage-Mackenzie (BLM) prescription [71] confirmed this result.
The γ * γ * cross section is obtained from Eq. (57) using the optical theorem. For transversely polarized photons the cross section is 3 while for the longitudinally polzrized gluons we get where we have also used 1 2 ln(a 1 a 2 ) = ln s Q1 Q2 along with the fact that F (ν) = F (−ν). The latter condition, while satisfied by F (ν) in Eq. (55), is also true in general, since it follows from the symmetry of the dipole-dipole scattering under the interchange of the dipoles. In Eqs. (58) and (59),Ĩ ρσ LO+NLO (q, ν) are given by (28) withᾱ µ → α s (Q 1 Q 2 ) and Re[F (ν)] is given by Eq. (55) [52]. Note thatĨ ρσ Eqs. (57), (58) and (59) are the main results of this Section, giving us the NLO forward scattering amplitude and the transverse (TT) and longitudinal (LL) cross sections for the γ * γ * scattering. Comparing Eqs. (58), (59) with (43), we make an important observation that the structure of the NLO γ * γ * cross section is the same as that at the LO with the addition of the running of the coupling and the corrections to the impact factors, intercept, and the "daughter" dipole-dipole scattering. The γ * γ * cross sections for different polarizations (e.g. TL or LT) can be constructed using Eq. (57) and the optical theorem.
are included into the γ * γ * → ρρ cross section, with Eq. (43) in [74] containing a power-of-energy form of the NLO evolution correction as a conjecture for how higher-order iterations of the NLO evolution corrections would come in. The power of energy in Eq. (43) of [74] agrees with the power of s/(Q 1 Q 2 ) in our Eq. (58), with the latter being an exact NLO BFKL eigenvalue resulting from construction of the NLO BFKL eigenfunctions in [46]. The γ * γ * → ρρ process considered in [74] is different (in the impact factors) from the γ * γ * → γ * γ * process that Eqs. (58) and (59) were derived for. Therefore, a complete side-by-side comparison of the two results is not possible at this point. Further work is needed to relate our approach in this paper and in [46] to that of [72][73][74].
D. Numerical result of the γ * γ * cross-section The NLO γ * γ * cross section we obtained in Eqs. (58) and (59) is an analytic result: it is a convolution of the two photon impact factors at NLO and the energy-dependent factor given by the exponentiation of the eigenvalues of the NLO BFKL equation. In this section we plot the result of the cross section as a function of rapidity Y = ln(s/Q 1 Q 2 ) and for several values of the virtualities of the two scattering photons (again omitting f e 2 f 2 ). We plot the γ * γ * cross section as a function of rapidity for virtual photons with momenta Q 1 = Q 2 = 1 GeV in Fig. 3. The left panel of Fig. 3 depicts the TT cross section from Eq. (58), while the right panel contains the LL cross section from Eq. (59). The one-loop coupling constant in our results (58) and (59) runs with the momenta of the virtual photons. In Fig. 3 we use Eq. (9) with N f = 3 and Λ QCD = 250 MeV; hence α s (1GeV 2 ) ≈ 0.50. Fitting the curves in Fig. 3 with e ∆ Y we make an estimate of the effective intercept; we get ∆ ≈ 0.78 for the TT cross section and ∆ ≈ 0.75 for the LL cross section. Compared to the LO BFKL intercept ∆ LO ≈ 4αs Nc π ln 2 ≈ 1.33 for α s ≈ 0.50 we conclude that the NLO BFKL evolution indeed tends to significanlty reduce the intercept compared to the LO result. However, as is clear from Fig. 3, the LO+NLO cross section keeps growing with rapidity even for rather large values of the coupling α s ≈ 0.50. This appears to be in qualitative agreement with the numerical solution of the NLO BFKL equation in [42]. Note also that the shape of our cross sections in Fig. 3 is similar to that for the γ * γ * → ρρ cross section from [74], though a detailed numerical comparison between the two can not be carried out due to different scattering processes considered. In Fig. 4 we plot the TT (left panel) and LL (right panel) cross sections for Q 1 = Q 2 = 5 GeV (dashed line) and Q 1 = Q 2 = 10 GeV (solid line). In Fig. 4 we use N f = 5 such that α s ((5GeV) 2 ) = 0.27 and α s ((10GeV) 2 ) = 0.22 respectively. The estimates of the intercepts are as follows: ∆ ≈ 0.34 for the TT cross section and ∆ ≈ 0.32 for the LL cross section for Q 1 = Q 2 = 5 GeV and ∆ ≈ 0.27 for both the TT and LL cross sections at Q 1 = Q 2 = 10 GeV. Comparing this to ∆ LO ≈ 0.73 for Q 1 = Q 2 = 5 GeV and ∆ LO ≈ 0.59 for Q 1 = Q 2 = 10 GeV we again see a reduction of the intercept at NLO, though again it remains positive corresponding to a cross section growing with Y .
We notice that the all the plotted cross sections are negative for small values of rapidity. Similar behavior has been observed in [74]. In fact, negativity takes place without BFKL evolution at Y = 0: we see that the NLO correction to the daughter dipole-dipole cross section (corresponding to the forward amplitude in Eq. (54)) is numerically large and negative. Comparing the curves in Figs. 3 and 4 we see that negativity tends to diminish with inclreasing Q/decreasing α s , which supports this interpretation. We have verified that at very large Q 1 , Q 2 the cross sections become positive-definite for all values of Y . We would like to point out that at small rapidities, Y ≈ 0, there are corrections to the dipole-dipole cross section which are inversely proportional to s and are given by the diagrams like that in Fig. 2b along with the subleading contributions to the diagram in Fig. 2a. Such corrections were neglected in deriving Eq. (54) [52]. However, such corrections are likely to become important near Y ≈ 0, modifying the dipole-dipole cross section and possibly making it positive again even for the range of Q 1 , Q 2 considered here. Hence the negativity of the cross sections in Figs. 3 and 4 at small Y is not necessarily a sign of the failure of the formalism, but rather an indication that energy-suppressed corrections need to be added in that region. Since the NLO cross sections in Eqs. (58) and (59) are likely to dominate over the energy-suppressed corrections at rapidities of the order Y ≈ 2 ÷ 3, we conclude that our NLO result gives a reliable prediction for the γ * γ * cross section only for sufficiently large values of Y , e.g. for Y > ∼ 2 in the kinematics of the plots in Figs. 3 and 4. Note also that at very high rapidity (Y = Y U ) the unitarity corrections due to nonlinear evolution will become important; hence the validity region of Eqs. (58) and (59) in rapidity is also limited from above, Note also that if the cross sections in Figs. 3 and 4 are replotted using Re [F (ν)] from Eq. (55) without the second term in parenthesis on its right-hand side (∼ 4/(4ν 2 + 1)), the existence of which was challenged in [69] as violating the principle of maximum transcedentality, they become positive-definite for all values of Y since in this case the NLO correction to the dipole-dipole cross section is also positive-definite. However, such cross sections are not monotonic functions of Y and tend to decrease at small Y before the growth picks up at larger Y .

A. Defining the BFKL Equation Beyond NLO
This Section is dedicated to constructing a solution to the NNLO BFKL equation in QCD. However, before we begin constructing the solution, it appears necessary to clarify what we mean by the BFKL equation beyond NLO. 4 Consider first a scattering of two quarkonia, say, resulting from γ * γ * decays as considered in the previous Section, mediated by the LO BFKL evolution. In terms of parameters α s and Y the corresponding cross section is proportional to with the last transition done in the power counting of Eq. (1) in which α s Y ∼ 1. The NLO cross section of Sec. III, with the NLO corrections included in the impact factors, evolution, and daughter dipole-dipole scattering is parametrically It is natural to assume that NNLO BFKL evolution, along with NNLO corrections to the impact factors and daughter dipole scattering would lead to the γ * γ * cross section of the form While this is the correct NNLO BFKL contribution to the γ * γ * cross section, there are in fact other corrections contributing to the γ * γ * cross section at the order-α 4 s in our power counting. The corrections result from multiple BFKL pomeron exchanges and from interactions mediated by the Bartels-Kwiecinski-Praszalowicz (BKP) multireggeon states [75,76]. These are schematically shown in Fig. 5 Single-splitting fan diagram or a single-loop diagram in Fig. 5 give the following parametric contribution to the onium-onium cross section [77,78] σ γ * γ * f an or loop ∼ α 2 s e const αs Y 2 ∼ α 4 s , which, in our counting, is of the same order as the NNLO BFKL (single-ladder) contribution to the cross section. The diagrams in Fig. 5 are, clearly, not included in the BFKL evolution equation. We conclude that the NNLO BFKL description of the dipole-dipole scattering receives order-1 corrections from the diagrams in Fig. 5, at least in terms of the α s and Y power counting (1). At the same time, one can define BFKL evolution with an arbitrary order kernel in α s by linearizing the arbitraryorder BK/JIMWLK evolution equations for scattering of a dipole on a dense target. In the large-N c limit, where BK equation for the dipole amplitude is valid to arbitrary order, mathematically it can always be linearized giving us the arbitrary-order BFKL equation. Beyond large-N c one can use the (first) equation in the Balitsky-JIMWLK hierarchy, the equation for the fundamental dipole: linearizing this equation one would still recover the arbitrary-order BFKL equation, now without the large-N c constraint. We see that mathematically one can always define BFKL evolution to an arbitrary order in α s as a linearization of the dipole evolution equation.
Note also that linearizing the higher-order evolution of the correlators of more than two Wilson lines, such as the fundamental quadrupole [79,80], and subtracting out the BFKL evolution pieces, allows one to also construct the BKP evolution equation [81] with the kernel calculated to any higher order in the coupling. Hence, in the s-channel language, the BFKL evolution is defined as the linearized evolution of a dipole, while the BKP evolution results from the linearized equation for correlators of higher number of Wilson lines (with BFKL exchanges subtracted out).
Indeed a physical question arises concerning the existence of a region where this linearization is justified, in the sense that the BFKL evolution obtained by BK/JIMWLK dipole evolution linearization dominates the scattering cross section. We can imagine two regimes in which an arbitrary-order BFKL evolution dominates over non-linear corrections.
Consider DIS on a dense target, in which saturation effects could be important. According to the conventional wisdom, the nonlinearities can be neglected for Q ≫ Q s (Y ), with Q s (Y ) the saturation scale describing the target for scattering with the total rapidity interval Y . (Note that saturation effects may still come in through the initial conditions to the linear BFKL evolution though [7,82,83].) For BFKL to be valid it is also important to make sure that the photon virtuality Q 2 is not large enough for the DGLAP evolution effects to become large. That is, we would want α s ln Q 2 /Q 2 s (Y ) ≪ 1. In the end we obtain the following range of Q in which the arbitrary-order BFKL evolution dominates the total DIS cross section For small enough α s this region could be large. In the dilute-dilute γ * γ * scattering case it appears to be less clear how to define a similar kinematic regime of BFKL dominance: in this case the large-N c limit comes in handy. The other regime in which the nonlinearities associated with multiple pomeron exchanges and other multi-reggeon states could be neglected is the very large-N c limit. When N c is very large, all interactions are suppressed, but the non-linear (or multi-reggeon) interactions are more suppressed than the linear BFKL interactions. Namely, a single-BFKL ladder exchange in the dipole-dipole scattering is parametrically of the order where we now keep track of the powers of α s and leading powers of Y and 1/N c . In the last step in Eq. (65) we have employed the 't Hooft large-N c limit. One can show that the multiple exchange diagrams, like those in Fig. 5, are suppressed by at least two powers of N c compared to the single BFKL ladder exchange (65). This is clear for the pomeron loop and fan diagrams (panels A and B) in Fig. 5, since they are explicitly non-planar, and, hence, N c -suppressed. For the BKP evolution in panel C the N c suppression becomes apparent after one subtracts out of it the single-and double-pomeron exchange contributions (see e.g. [81]). In the end we see that and these corrections can be neglected when compared to the single BFKL ladder exchange in the large (but finite) N c limit. This is why in the large-N c limit BFKL equation is dominant up to any order in the coupling α s . This conclusion is confirmed by the fact that it makes sense to define BFKL equation even at strong coupling, as shown in [84,85] using the anti-de Sitter space/conformal field theory (AdS/CFT) correspondence [86,87]. (Note that BK evolution, while being large-N c in the projectile wave function, does resum subleading-N c corrections in the multiple interactions with the target: while being 1/N 2 c -suppressed, such corrections are enhanced by powers of A 1/3 for a nuclear target with atomic number A, or, equivalently, by powers of color charge density ρ of partons in a proton or a nucleus.) To summarize, we have argued that it is possible mathematically to define the BFKL equation at any order of the perturbation theory as the linearization of the dipole evolution equation. Moreover, we have identified two physical situations in which the BFKL exchange dominates over other interactions. This clarifies what we mean by solving the BFKL equation at NNLO, as we will do in the remainder of this Section.

B. Eigenfunctions and Eigenvalues of the NNLO BFKL Equation
In this Section we derive the eigenfunction of the NNLO BFKL equation in the n = 0 azimuthally-symmetric case. Generalization to n = 0 case can be easily accomplished along the lines of Sec. II.
Similarly to the NLO eigenfunction case we obtained in our previous paper [46], we look for NNLO BFKL eigenfunctions by adding perturbative corrections to the LO eigenfunction up to NNLO in the following way: The NLO correction to the BFKL eigenfunctions was found in [46] yielding (cf. Eq. (13)) Our goal is to construct the NNLO correction f to order-α 3 µ . (The full BFKL kernel is K(k, q) = K LO+NLO+NNLO (k, q) + O(α 4 µ ).) Indeed the full NNLO BFKL kernel is not known at the time of writing this paper. However, we will not need to know the NNLO BFKL kernel in order to find its eigenfunctions. In [46] we have constructed the NLO BFKL kernel's eigenfunctions employing only the projection of the NLO BFKL kernel onto the LO BFKL eigenfunctions (powers of momentum). Similarly here we will only need a projection of the NNLO BFKL kernel onto the LO BFKL eigenfunctions. While this projection is also not known in general, its structure could be inferred from that of the NLO kernel projection (7): the NNLO projection should consist of an order-α 3 µ conformal part along with the k-dependent one-and two-loop running coupling corrections to the LO+NLO projection (7). We thus write (cf. Eq. (63) in [46]) Here the two-loop QCD beta-function is given by An important point in the determination of the NLO BFKL eigenfunctions in [46] was to show that the corresponding eigenvalues were real since the NLO kernel acts as a Hermitean operator on its eigenfunctions. For this reason, although we do not need to know exactly the expression of the NNLO BFKL kernel, we have to split the conformal part of its projection on to the LO eigenfunctions into the imaginary (odd under the ν → −ν) and real (even under ν → −ν) contributions. Thus, we define the real part of the conformal contribution as χ 2 (ν) and the imaginary part one as i δ 2 (ν) in Eq. (70). Let us stress again that at this point neither χ 2 (ν) nor δ 2 (ν) are known: however, at the end of this calculation we will obtain an explicit expression for δ 2 (ν) (see Eq. (96) below).
Let us consider the action of the NNLO BFKL kernel on the NNLO eigenfunction (67). We require that the eigenfunction condition is satisfied, with ∆(ν) denoting the BFKL eigenvalue. The LO and NLO contributions to the eigenvalue are known [33,46] (see (16)). We parametrize the NNLO contribution to the eigenvalue in terms of another unknown function c (2) (ν) such that Our task is to substitute Eqs. (67) and (73) into Eq. (72) and find f (2) ν (k) and c (2) (ν) satisfying the resulting equation at order-ᾱ 3 µ . We begin by evaluating the left-hand-side of Eq. (72) at order-ᾱ 3 µ : using Eqs. (67) and (69) we write The first term on the right-hand-side of (74) can be read off from Eq. (70): To find the second term on the right-hand-side of (74) we write, using Eq. (68), The projection of K NLO onto q −1+2 i ν can be read off from the order-ᾱ 2 µ terms on the right side of Eq. (70). Using that in Eq. (76), after some algebra one arrives at Henceforth, for brevity, we will use χ 0 ≡ χ 0 (ν), To proceed we now need an ansatz for f (2) ν (k): similarly to the NLO case [46], we expand it in the powers of ln k 2 /µ 2 , In the NLO case we truncated the series at n = 2. What determines the order at which we should truncate the series is the minimum number of terms needed in order to make the terms proportional to powers of ln k 2 /µ 2 in equation (78) disappear. A quick analysis of Eq. (78) shows that in the NNLO case one has to truncate the series at n = 4. We thus write 1 (ν) ln The action of the LO BFKL kernel on the NNLO part of the eigenfunction can be found using the standard replacement of logarithms by derivatives, Substituting Eqs. (80) and (81) into Eq. (78) we obtain (with c By equating the powers of ln(k 2 /µ 2 ) in Eq. (82) we are able to determine c 3 and c 3 (ν) = iβ 2 2 − 5 12 and also find an equation which relates the coefficient c 1 to the function c (2) (ν): To further fix c 1 and c (2) we need to make sure the eigenfunctions H 1 2 +i ν (k) form a complete set of functions. As we have already shown in [46] the functions H 1 2 +i ν (k) have to satisfy the completeness relation order-by-order in the perturbation theory. This means that since the LO eigenfunctions already satisfy this condition, then all the higher-order inᾱ µ corrections to the eigenfunctions have to cancel on the left-hand-side of Eq. (85), possibly giving us an extra constraint on c 1 and c (2) . At order-ᾱ µ , the eigenfunctions H 1 2 +i ν (k) from (67) already satisfy the completeness relation provided that f (1) ν (k) is given by (68). At order-ᾱ 2 µ the completeness relation leads to the following constraint Plugging in f To see this one has to again use the relation ln(k 2 /k ′2 ) (k 2 /k ′2 ) i ν = −i ∂ ν (k 2 /k ′2 ) i ν together with integration by parts to eliminate terms containing powers of ln(k 2 /k ′2 ), such that only terms containing ln(kk ′ /µ 2 ) are left. Using this technique one then has to require that the coefficient of each power ln n µ 2 (with n = 1, 2, 3, 4) is zero independently; it turns out that all the terms proportional to ln 4 µ 2 , ln 3 µ 2 , ln 2 µ 2 are identically zero, while the terms proportional to ln µ 2 give the condition (87) for Re[c On the other hand the terms without ln µ 2 in (86) result in the following condition for Re[c (2) 0 (ν)] and Im[c where χ (5) 0 = ∂ 5 ν χ 0 (ν). We can now use the result (88) for Re[c (2) 1 (ν)] in (84) to determine Im[c (2) and also obtain a relation between Re[c (2) ](ν) and Im[c (91) (2) ν (k) and c (2) (ν), then introduced an ansatz for f (2) ν (k) which effectively expressed all the unknowns in terms of a total of six complex-valued functions: c 4 (ν) and c (2) (ν). Then, by requiring that the functions H 1 2 +i ν (k) are eigenfunctions of the NNLO BFKL kernel and that they satisfy the completeness relation (85), we were able to determine c 4 (ν) given by Eq. (83), along with Re[c (2) 1 (ν)] given in Eq. (88) and Im[c (2) (ν)] given in Eq. (90). We have also obtained equations (91) and (89) relating Re[c (2) (ν)] to Im[c (2) 1 (ν)] and Im[c 0 (ν). This freedom of choosing any c (2) 0 (ν) is analogous to the same freedom observed in constructing the NLO BFKL eigenfunctions in [46]: as we argued in [46], this is a consequence of the freedom of choosing any overall phase for the eigenfunctions along with the freedom to reparametrize the ν variable (by shifting it by an O(ᾱ 2 µ ) correction at NNLO). We can now write down the expression for the NNLO corrections f (2) ν (k) to the BFKL kernel eigenfunction (we use, as usual, the short hand notation where 1 ] ln The NNLO eigenfunctions are then given by where f (1) ν (k) is given in Eq. (68) and f (2) ν (k) is given in Eq. (92). Before we proceed, let us note that the eigenfunctions H 1 2 +i ν (k) satisfy the orthogonality expression as well, as can be verified explicitly to orderᾱ 2 µ . The BFKL kernel eigenvalues up to NNLO corrections are where we have used Eq. (91) for Re[c (2) (ν)] and Eq. (90) for Im[c (2) (ν)]. The BFKL kernel has to be hermitean at any order in the coupling constant: therefore, the BFKL kernel eigenvalues have to be real. Requiring the eigenvalues (95) to be real we fix the imaginary part i δ 2 (ν) of the conformal piece of the projection of the NNLO BFKL kernel onto the LO eigenfunctions. We thus predict that Note that this result is in complete agreement with our NNLO ansatz in [46] (see Eq. (B6) there). Thus, the NNLO BFKL kernel eigenvalues reduce to ∆(ν) =ᾱ µ χ 0 (ν) +ᾱ 2 µ χ 1 (ν) +ᾱ 3 µ χ 2 (ν) + Im[c At this point we have NNLO BFKL eigenfunctions (92) and the eigenvalues (97) for any function c (2) 0 (ν). One could use the phase-choice and ν-reparametrization freedoms to select a convenient value of c (2) 0 (ν). Alternatively one can argue that, since the functions H 1 2 +i ν (k) form a complete and orthogonal set of eigenfunctions of NNLO BFKL kernel for any c (2) 0 (ν), a particularly simple choice is to put Im[c (2) 0 (ν)] = 0, Im[c with Re[c (2) 0 (ν)] easily found from Eq. (89). We then obtain our final expression for the NNLO BFKL eigenvalues and eigenfunctions The NNLO BFKL eigenfunctions are given by Eq. (93) with f Notice that the only unknown function in (99) is χ 2 (ν) which can be found only after performing the calculation of the NNLO BFKL kernel (with NNLO BFKL evolution defined in Sec. IV A). Steps in this direction were done in [88,89].
Another important observation is that to determine the structure of higher-order corrections to the BFKL eigenfunctions and eigenvalues, one only needs to know the running coupling part of the projection of the corresponding higher-order BFKL kernel onto the LO eigenfunctions (see, e.g., Eq. (70)), which can always be obtained as runningcoupling corrections to the lower-order projections. Hence our procedure can be used to construct NNNLO and higher-order BFKL eigenfunctions, and is limited only by our knowledge of the QCD beta-function.

C. General Form of the NNLO BFKL Equation Solution
In [46] we have obtained a solution of the NNLO BFKL equation for the Green function (2) (provided that χ 2 (ν) is explicitly calculated), in an indirect way: we assumed that the NNLO solution could be obtained by a simple generalization of the NLO expression for the BFKL Green function to the NNLO order. It turned out that such an ansatz did not solve the NNLO BFKL equation and had to be augmented to solve NNLO BFKL with a particular choice of δ 2 (ν), incidentally the same as in Eq. (96). Our goal now is to construct the NNLO expression for the BFKL Green function and to verify the ansatz proposed in [46] (see Eq. (B7) there).
The BFKL Green function can be written in terms of eigenfunctions as with H 1 2 +i ν (k) given by Eqs. (68), (100), (93), and ∆(ν) given in Eq. (99) at NNLO. Substituting those results in Eq. (101), employing the same trick of turning each ln(k 2 /k ′2 ) into −i ∂ ν acting on the powers of transverse momentum and integrating by parts, after a lengthy calculation one can rewrite Eq. (101) as with ∆(ν) given in Eq. (99). Using the fact that, at NNLO accuracy, along with we can write the NNLO BFKL solution as This result exactly coincides with, and, therefore, proves the ansatz suggested in Eq. (B7) of our previous paper [46]. Note that the first term in the square brackets of Eq. (105) taken in the saddle-point approximation of LO BFKL equation (that is, at ν = 0), gives exactly the ∼ α 5 s Y 3 term in the exponent found in [36]. With the accuracy of our NNLO approximation we can not distinguish the O(α 2 s )-terms in the exponent from those in the prefactor: it is, therefore, possible that the first term (along with the other terms) in the square brackets of Eq. (105) would exponentiate when high-order corrections are included. An exponential form of the first term in the square brackets of Eq. (105) without the saddle point approximation was obtained in [47].

V. CONCLUSIONS
There are two main results presented in this paper. The first one is the analytic NLO γ * γ * cross section given in Eqs. (58) and (59)  The structure of the NLO γ * γ * cross sections (58) and (59) is the same as the one at LO: the two impact factors, the dipole-dipole cross section and the energy-dependent exponential factor each receive corrections at NLO; in addition each coupling constant is replaced by a running coupling. The factorization of the scattering amplitude is represented in Fig. 1 and is likely to persist at higher orders (in the linear approximation). The corresponding expression (57) for the forward γ * γ * scattering amplitude, involving projections of the impact factors on the NLO eigenfunctions, is also likely to be true to all orders.
The γ * γ * cross sections (58) and (59) plotted in Figs. 3 and 4 are negative at small values of rapidity Y . This issue may be resolved by adding the energy-suppressed contributions to the cross sections, which become important at low energy (rapidity Y ) making the net cross section positive for all Y .
An ansatz for the structure of the solution of the NNLO BFKL equation was constructed [46], but it was obtained in an indirect way by attempting to guess the NNLO BFKL solution using an analogy to the NLO solution. In addition the ansatz from [46] relied on a particular expression for δ 2 (ν). In this work we have obtained the solution of the NNLO BFKL equation by explicitly constructing the NNLO eigenfunctions in (93). This procedure yielded the NNLO BFKL Green function which confirmed the ansatz from our previous paper [46]. Finally, Fourier-transforming Eq. (A7) into transverse momentum space with the help of Eq. (26) yields For the use in the main text, let us define and observe that f (ν) f (−ν) = 4 ν 2 (4 ν 2 + 1) 2 .
With the help of definition (A9) we rewrite Eq. (A8) as Appendix B: NLO γ * γ * cross section In this section we provide details of the calculation leading to Eqs. (54), (56) and (57). In the derivation we will make use of the solution of the NLO BFKL equation which, replacing the logarithm ln k with −i∂ ν acting on k iν and integrating by parts, can be written as Eq. (B2) is a result we obtained in Ref. [46]. We proceed with Eq. (54). Using Eqs. (52) and (46) we write Employing Eqs. (A11) and (45) we write Using the result of [62] (or the calculation of [52] modified for QCD) we write (cf. Eq.