Solution to the Balitsky-Kovchegov equation with the collinearly improved kernel including impact-parameter dependence

The solution to the impact-parameter dependent Balitsky-Kovchegov equation with the collinearly improved kernel is studied in detail. The solution does not present the phenomenon of Coulomb tails at large impact parameters that have affected previous studies. The origin of this behaviour is explored numerically. It is found to be linked to the fact that this kernel suppresses large daughter dipoles. Solutions based on a physics motivated form of the initial condition are used to compute predictions for structure functions of the proton and the exclusive photo- and electroproduction of vector mesons. A reasonable agreement is found when comparing to HERA and LHC data.


I. INTRODUCTION
Evolution equations are powerful tools to study the high-energy, equivalently, smallx limit of quantum chromodynamics (QCD) [1][2][3]. The availability of quality data from HERA [4] and the LHC [5] as well as the need for reliable phenomenology for the proposal of new electron-ion facilities [6,7] have given an extra impulse to the development of these tools.
Soon after its introduction, the kernel of the leading order BK equation was modified to include corrections that take into account the running of the coupling constant [20][21][22]. The resulting equation, referred to as rcBK below, when combined with appropriate initial conditions -embodying non-perturbative properties of the hadronic targets -and disregarding the impact-parameter dependence, produces solutions that have been successfully used to described a wide variety of phenomena. In particular, the structure function data of the proton as measured at HERA was successfully described [23][24][25][26]. A few other applications of these solutions are, for example, gluon production in heavy-ion collisions [27], single particle [28] and J/ψ production in pp and pA collisions [29] , di-hadron correlations in p-Pb interactions [30] and even the flux of atmospheric neutrinos [31,32].
As already mentioned, these comparisons of rcBK-based predictions to data disregarded the impact-parameter dependence of the dipole amplitude. The reason is that earlier studies of solutions including the impact parameter found that the amplitude developed a power-like dependence on b ≡ | b|, the so-called Coulomb tails, which generate an unphysical growth of the cross section [33]. Nonetheless attempts were made to modify the kernel to solve this problem, for example, by adding an ad-hoc cut-off for large sizes of the daughter dipoles [34].
The solutions found had no more Coulomb tails, but needed an extra, so-called soft, contribution to be able to describe HERA data on structure functions [35]. (A similar conclusion also holds for the solutions of the impact-parameter dependent JIMWLK equation [36].) Nonetheless, this approach did a good job when confronted with HERA data on exclusive vector meson production [37].
Recently, the kernel of the leading order equation has been improved by including the resummation of all double collinear logarithms [38] as well as two classes of single logarithmic corrections [39]. Using this kernel and disregarding the dependence on the impact parameter, it was also possible to obtain a good description of HERA data on the structure function of the proton. Finally, in the rapid communication [40], we have demonstrated that solutions of the BK equation with the collinearly improved kernel and an appropriate initial condition describe correctly the HERA data on structure functions and the t dependence of the exclusive photoproduction of J/ψ at one energy without the need of any additional ad-hoc parameter or correction.
In this contribution the studies reported in [40] are extended to discuss in depth the behaviour of the collinearly improved kernel and of the solutions of the corresponding BK equation, comparing them to the rcBK case. In addition, more details on the comparison to HERA structure function data are presented, and comparison of our predictions to relevant HERA and LHC data on exclusive vector meson photo-and electro-production is provided.
In all cases, the agreement between model and measurements is satisfactory.
The rest of this contribution is organised as follows: In Sec. II the formalism used throughout this work is reviewed. In Sec. III the technical details to solve the collinearly improved impact-parameter dependent BK equation are addressed. In Sec. IV the origin of the suppression at large impact parameters is discussed, the behaviour of the solution is contrasted with solutions of the rcBK case, and the shape of the amplitude is shown at different values of rapidity, dipole size and impact parameter. In Sec. V and Sec. VI our predictions are confronted with structure function data measured at HERA, and to data for cross sections of exclusive photo-and electroproduction of φ, J/ψ, ψ(2S), and Υ(1S) vector mesons measured both at HERA and at the LHC, respectively. Sec. VII contains a brief summary of our findings and presents our conclusions.

A. The Balitsky-Kovchegov equation
The BK evolution equation reads [20,21] where r ≡ | r |, r 1 ≡ | r 1 |, and r 2 ≡ | r 2 | ≡ | r − r 1 | are the sizes of the original dipole and of the two daughter dipoles, respectively. Note that these are 2-dimensional vectors in the same plane as the impact parameter. The magnitudes of the corresponding impact parameters In this work, the solution to the BK equation is obtained under the assumption that the scattering amplitude N( r, b, Y ) depends solely on the sizes of the dipoles and of the impact parameter vectors. In practice, this means to solve the equation subjected to the condition that the angle between r and b is fixed. We chose to fix this angle at zero, meaning that these vectors are parallel.

B. Kernels of the Balitsky-Kovchegov equation
Several functional forms for the kernel of the BK equation have been proposed. The ones that are mentioned in this work are presented in the following.
The leading order kernel is given by where the non-running coupling, α nr s , is fixed to a constant value. The running coupling kernel K rc (r, r 1 , r 2 ) reads [20] K rc (r, r 1 , where N c is the number of colours and α s is the running coupling, which is further discussed in Sec. II C.
The running coupling kernel with a cutoff to tame the Coulomb tails generated by the evolution in the impact parameter is given by [35] K bdep rc (r, r 1 , r 2 ) = K rc (r, r 1 , where Θ is the Heaviside function and m a parameter to limit the size of daughter dipoles. Finally, the collinearly improved kernel is [39] K ci (r, r 1 , r 2 ) = α s 2π where J 1 is the Bessel function (the inclusion of the Bessel function into the BK kernel has been previously discussed in [41]), the anomalous dimension is A 1 = 11/12, and The sign factor in the exponent ±α s A 1 takes the value of the plus sign when r 2 < min(r 2 1 , r 2 2 ) and the negative sign otherwise. For the running coupling the smallest dipole prescription is used throughout the computation according to where r min = min(r 1 , r 2 , r). This prescription was compared to other prescriptions in [39], where it was found to work adequately in this context. This prescription has also been suggested as the natural option for the BK equation at next-to-leading order [42].

C. Treatment of the coupling constant
In this work the running coupling is computed in the variable-number-of-flavours scheme, implemented according to where n f corresponds to the number of flavours that are active, C 2 is an infrared regulator that takes into account the approximations made for the computation of the Fourier transform into the position space and is usually fit to data [24]. The variable β 0,n f is the leading order coefficient of the QCD beta-series and is given by relation The value of the QCD scale parameter Λ 2 n f depends on the number of active flavours. When heavier quarks are active (charm and beauty quarks), its value is obtained from the rela- This recursive relation needs to be fixed at one point and for this the usual choice is to take the value of the running coupling at the scale of the mass of the Z 0 boson. In this way, Λ 5 is set with the use of the experimentally measured value of α s (M Z ) = 0.1196 ± 0.0017, where the Z 0 mass is M Z = 91.18 GeV/c 2 [43]. The number of active flavours is set depending on the transverse size of the mother dipole. The condition that governs this relates the mass of the heaviest quark considered to the values of the dipole size r. This condition can be expressed as Since all dipole sizes are accounted for in the BK evolution equation, there is a need to freeze the coupling at a set value after a certain dipole size is reached [24]. In this work, the coupling is frozen at α sat s = 1 as in [38].
The value of the parameter C affects the description of data by modifying the speed of the evolution and effectively changes the slope of the structure function. The higher value of this parameter the more the running of the coupling is suppressed and, consequently, the slope in the structure function F 2 is less steep. Figure 1 compares the running of α s for two values of C: the one used here, C = 9, and the one used in [38], C = 2.586. The value C = 9 was set heuristically and since the solutions reproduce correctly the data, as shown below, it has not been further optimised. with C = 2.586 (red) and C = 9 (blue).

A. Initial condition
The initial condition, already introduced in [40], depends on the impact parameter; it is suppressed in the regions where r or b reaches large values, in order to respect the geometric nature of the dipole-proton interaction. The shape of its functional form is a combination of the expected behaviour in r, which is obtained from the Golec-Biernat Wüsthoff (GBW) model [44], and the impact-parameter dependence, which uses a Gaussian distribution to reflect the expected profile of the proton. Such an approach has been used in similar forms in the past; e.g., in [45][46][47][48][49]. The main new ingredient with respect to the initial condition used in the previous studies [34,35,37] is the explicit separation of the contribution from the individual quark and anti-quark forming the dipole. The initial condition is given by where b q i are the impact parameters of the quark and anti-quark forming the dipole and As a first attempt, the angle between r and b was fixed as shown schematically in Fig. 2.
As the results obtained with this initial condition are satisfactory, no further optimisation has been considered.
The parameters appearing in this initial condition, Q 2 s and B G , have a clear physical interpretation as the saturation scale and as the variance of the Gaussian distribution of the target in impact parameter, respectively. The value of the Q 2 s parameter is chosen to be 0.496 GeV 2 , such that the F 2 (x, Q 2 ) data are correctly described at the rapidity of the initial condition. The relation between x and rapidity is Y = ln x 0 /x , where x 0 = 0.008. The parameter B G is set to 3.2258 GeV −2 in order to describe the data for exclusive photoproduction of J/ψ off protons at a photon-proton centre-of-mass energy W γp = 100 GeV, where, as customary x = (M 2 + Q 2 )/(W 2 γp + Q 2 ) is used; here, M represents the mass of the vector meson.

B. Setup for the numerical solution to the equation
The BK evolution equation does not have an analytic solution and therefore has to be solved numerically. The procedure used by us in [26,50] was extended to the case of the impact-parameter dependent BK equation [40] and the solution is evolved in rapidity with a step of ∆Y = 0.01. Since the transverse dipole vectors are related as r = r 1 + r 2 , by fixing the values of r and r 1 to the pre-defined grid, the values of r 2 are in general off the grid. Whenever this happens, linear interpolation in the log 10 space is used to get the desired value of N(r 2 , b 2 , Y ).
A similar approach is used for obtaining the value of the scattering amplitude whenever the value of b 1 or b 2 is off the grid.
The values of b 1 and b 2 are then computed from the relations b 1 = b+ r 2 /2 and b 2 = b− r 1 /2 assuming a fixed angle between r and b. As mentioned above, this angle is set to zero for the results presented below.
The solution to the BK equation has been implemented independently using C++ and the Intel Fortran Compiler. Both implementations have similar performance, with the Fortran version being slightly faster. In a standard personal computer, the program performs the evolution of the dipole amplitude in one unit of rapidity, that is 100 steps for the settings described above, in a bit less than one hour for one set of parameters.
To test the numerical stability of the selection of the grid, the setup was modified and the scattering amplitude was compared at Y = 3, r = 1/GeV and all values of b. We have changed the step in rapidity from 0.01 to 0.02, the number of steps in r and b per order of magnitude from 25 to 15 and the size of the grid in the polar angle from 21 to 16 and 31 points. Except for the change to 16 points in the grid for polar angles, all other changes produced a difference below the per-mil level. The use of the spare grid in polar angle produced changes almost at one percent level. In summary, with the chosen settings a numerical precision at the percent level, or even below it, is expected.

IV. THE SOLUTION TO THE BK EQUATION
A. Behaviour of the collinearly improved kernel As was shown in [40], the solutions to the BK equation do not exhibit Coulomb tails when using the collinearly improved kernel. This behaviour is related to the suppression of this kernel for large values of the size of the daughter dipoles. As an illustration, Fig. 3 shows the ratio of the collinearly improved kernel, see Eq. (4), to the running-coupling kernel, see Eq. (6). (The parameter C for the running coupling in this kernel was chosen to be C = 9 just as in the collinearly improved kernel for the sake of a valid comparison.) The ratio is computed at r = 1 GeV −1 and θ rr 1 = π/2. Other values produce a similar picture. The figure shows that for large sizes of the daughter dipole the collinearly improved kernel is orders of magnitude smaller than the running-coupling one.
To follow up in more detail the origin of this behaviour the kernels are divided into three parts. For the collinearly improved kernel, they are The first term, K 1 ci , is present already at the leading order if one considers a fixed value of the running coupling, K 2 ci takes into account the contribution from the single collinear logarithms, and K 3 ci resumms double collinear logarithms to all orders. The entire collinearly improved kernel is then given by the multiplication of all these factors as For the running coupling BK kernel, the separation in three parts is as follows: whereas the running coupling kernel is then given by the addition of these constituent terms as The contribution of the three terms is shown in Fig. 4 at r = 1 GeV −1 and θ rr 1 = π/2 for each of the two kernels. The fact that the three terms are added in K rc , but multiplied in K ci explains numerically the suppression. Even though, the first term is essentially the same for both kernels, the additive character of K rc makes it deviate from the collinearly improved kernel at large r 1 values as shown in Fig. 4. There, we can see that even though the kernels are comparable in the low-r 1 region, at large r 1 values, the K 2 rc and K 3 rc terms become dominant, whereas in the collinearly improved kernel, the K 1 ci term suppresses the total value. contribute to the total integral in Eq. (2). This is true because a large impact parameter means that the probing dipole is far away from the target proton and the amplitude is therefore (at the initial condition) exponentially suppressed. Only dipoles with r 1 (r 2 ) ∼ 2b can be oriented so that their impact parameters b 1 (b 2 ) are small, such that they contribute to the evolution. But, since K ci is suppressed in this region, the integral will be suppressed as well and the scattering amplitude will not grow fast at large b.
This can be numerically studied by computing the contribution to the evolution of the three terms in the collinearly improved kernel. Figure 5 shows the scattering amplitude after evolution to Y = 3 using each time a kernel formed with different constituents. It is clearly seen that the impact parameter profile is mostly influenced by the inclusion of the K 3 ci term with the Bessel functions. This term originates from resumming double collinear logarithms. Note that also the term K 2 ci , resumming single collinear logarithms, suppresses the large b region. initial condition is lifted with evolution. Similar behaviour was observed in previous studies, e.g. [34]. Nonetheless, in the case of the collinearly improved kernel the growth at the largest dipole sizes is not as fast and a shoulder appears, after which the amplitude is again suppressed. The behaviour as a function of impact parameter has been discussed above; the profile impact parameter grows, but the development of Coulomb tails is suppressed.
Recently, a similar finding has been reported for the case of NLO BFKL equations at large impact parameters [52].
Finally, Fig. 8 shows N(r, Y ), defined as for different dipole sizes and for two kernels, the running coupling and the collinearly im- proved. For small dipoles the difference is larger and it grows with rapidity. At larger dipole sizes the difference between both kernels is smaller. Note that for the comparisons to data discussed below, the main numerical contribution comes from the region of relatively large dipoles. For the case of the structure function the main contribution for virtualities between 1 and 10 GeV 2 comes from dipoles of sizes on the range around 0.1/GeV to 10/GeV, see e.g. the lower panel of Fig. 4 in [26].
Another interesting observation is that N(r, Y ) is related to the σ 0 parameter introduced in studies based on the rcBK equation without impact-parameter dependence. Basically, σ 0 corresponds to the scale of N(r, Y ). Standard values found for this parameter are a few tens of mb, see e.g. Table 1 in [25]. notation of [44] are used: and where m q i is the mass of the considered quark, which is set to 100 MeV/c 2 for light quarks and 1.3 GeV/c 2 for charm quark and 4.5 GeV/c 2 for bottom quark. Note that the computed structure function does not depend strongly on the value of the mass of the light quarks (as was reported in [39]); this has been checked by also using m u,d,s = 10 MeV/c 2 , which did not influence the description of data. The total wave function then is According to the optical theorem, one can link the dipole-target cross section to the scattering amplitude by dσ qq ( r, x) Furthermore, it is usual to shift the value of the x at which the structure function and reduced cross section are computed according to the photoproduction kinematic shift [44] x = x 1 + 4m 2 Using these ingredients, the relation for the computation of the structure function in the dipole model framework is and the reduced cross section is computed as where y = Q 2 /(sx) is the inelasticity of the process, s is the squared of the centre-of-mass energy of the collision and F L (x, Q 2 ) is given by the relation

B. Comparison to HERA data
The predictive power of this model is evaluated by confronting it with data from HERA on the F 2 (x, Q 2 ) structure function [53] in Fig. 9. A closer look is given in Fig. 10 for two values of the photon virtuality. To quantify the level of agreement between data and model, Fig. 11 presents the percentage pulls associated with the structure function, which are given by and by D % , which denotes the average of the corresponding values of d % . Finally, for completeness Fig. 12 shows the comparison of the model and data for the charm component of the proton structure function measured at HERA [53].
Overall, the agreement between prediction and data is within a few percent over most of the phase space. For our purposes this level of agreements is satisfactory. First, the equation we are using does not include the full angular dependence. Second, we have not needed to add any ad-hoc component to describe data and the values of the parameters are reasonable from the point of view of the physics that is being probed. Furthermore, note that the BK equation that we are using is definitely not the last word on the subject. The full equation at NLO has already been computed [42], and a large effort is being done to use it for phenomenology [51,[54][55][56][57]. There are also recent developments regarding the most adequate variable to evolve the scattering amplitude [58].

A. Exclusive cross section in the colour-dipole approach
Similarly to the DIS process described in the previous section, the diffractive production of a vector meson as a result of the interaction of a virtual photon with the proton can be treated within the colour-dipole approach. In this formalism, the exclusive cross section to produce a vector meson V is given by   where A T,L is the scattering amplitude of the process. It is given as a convolution of the overlap of photon-meson wave functions with the dipole cross section given in Eq. (30) (for a detailed derivation see e.g., [59]) and takes the following form where the subscripts T , L denote the contribution from the virtual photon with transverse, respectively longitudinal, polarisation, Ψ γ * is the wave function of a virtual photon which fluctuates into a dipole, Ψ V represents the wave function of the vector meson, and ∆ 2 ≡ −t, the square of the four momentum transferred in the proton vertex. Under the assumption of large photon-proton centre-of-mass energies W γp , where M is the mass of the given vector meson.
The wave functions of a vector meson are modelled under the assumption that the vector meson is predominantly a qq pair with the same polarisation and the spin structure as the original photon. The overlap of the photon-meson wave functions in Eq. (37) is given as and withê f being the effective charge of the given vector meson, ǫ defined by Eq. (28), and the parameter δ is a switch to include (δ = 1) or exclude (δ = 0) the non-local term in the longitudinal contribution. The scalar part φ T,L of the wave function is, in general, model dependent. For our studies, we use the boosted Gaussian model [60][61][62] in which the δ parameter is fixed to one. The values of the parameters for the wave functions of all vector mesons are fixed according to Table I in [63]. for details see [59]. The other correction takes into account that there are two values of x involved in the interaction of the dipole with the proton and one should therefore use the offdiagonal gluon distribution for vector meson production. This effect can be accounted for by multiplying the scattering amplitude by a factor R T,L g , called the skewedness correction [64].

B. Comparison to data
Using the model described in this paper, the cross sections for exclusive photo-and electroproduction of φ, J/ψ, ψ(2S), and Υ(1S) vector mesons are presented at different virtualities of the exchanged photon and they are compared to available experimental data.
The presented results are calculated at the scales which allow perturbative treatment of the specific parts of the model.
In Fig. 13   and ZEUS [66] for the |t| dependence (left) and the W γp dependence (right) of the exclusive electroproduction cross section of the φ meson. give very good agreement with the data at energies W γp = 55 GeV and W γp = 100 GeV.
The result for W γp = 78 GeV is slightly underestimated at low values of |t|, however one can notice the very small difference in the measured data with respect to the result for W γp = 100 GeV. Since the value of W γp from the experimental data is a mean value estimated from a measured energy range, the result of the model can be considered satisfactory.   [67,68] and LHC data from ALICE [69,70] for the W γp dependence of the exclusive photo-and electroproduction cross section of the J/ψ meson (left) and with HERA data from H1 [71] and ZEUS [72] for the W γp dependence of the exclusive photo-and electroproduction cross section J/ψ/ψ(2S) ratio (right).
The same comparison for the electroproduction at three different values of Q 2 can be seen in the right panel of the same figure. Although our predictions do not describe all the data points, we conclude the agreement between the data and the model to be qualitatively good.
The same conclusion applies to the comparison of the model predictions with the measured W γp dependence of the exclusive differential photo-and electroproduction cross sections at several fixed values of |t| presented in Fig. 15. The agreement of the predictions with the data is very good at low values of W γp , however at larger values (∼ 10 2 GeV), the predictions are underestimated when compared to experimental photoproduction data. We have also obtained total cross section for the J/ψ production which is presented in the left panel of the Fig. 16. The predictions for the electroproduction at three different values of Q 2 give a very good description of the available data. The result for photoproduction gives a good agreement with the data at low values of W γp , however at high energies the result is again underestimated when compared to data.
Also, the exclusive cross section of the ψ(2S) meson was calculated within the model.
The experimental data are not available for the total cross sections, but only for a ratio of the ψ(2S) to J/ψ cross sections, the predictions for these ratios for photoproduction and electroproduction at Q 2 = 16 GeV 2 are calculated and compared to data from H1 [71], and ZEUS [72], respectively, in the right panel of the Fig. 16. The description of the data is not very good, yet the large uncertainties of the experimental data do not allow us to make any final conclusions in this case.
To complete the set of the predictions based on the BK equation, the exclusive photoproduction of the Υ(1S) meson is presented in Fig. 17. The prediction is compared with experimental data obtained at HERA by H1 [73] and ZEUS [74] experiments. It is also compared with the two latest measurements -in proton-proton collisions at √ s = 7 TeV and √ s = 8 TeV by LHCb [75], and in proton-lead collisions at √ s = 5.02 TeV by the CMS experiment [76]. The description of the data is good, although the large uncertainties prevent us from making any strong conclusions regarding the agreement of the predictions with the data.

VII. CONCLUSIONS
The solution of the Balitsky-Kovchegov equation with the collinearly improved kernel and including the impact-parameter dependence has been obtained numerically. This solution does not show the so-called Coulomb tails that have appeared in previous attempts to include the impact-parameter dependence. We have shown that the suppression at large values of the impact parameter is due to the suppression of contributions from daughter dipoles of large sizes in the terms of the collinearly improved kernel that deal with the resummation  [73] and ZEUS [74], and LHC data from LHCb [75] and CMS [76] for the W γp dependence of the exclusive photoproduction cross section of the Υ(1S) meson.
of double and single collinear logarithms.
The solutions based on a physics-inspired initial condition have been confronted with HERA and LHC data of the structure function of the proton measured in deep-inelastic scattering and of exclusive vector meson photo-and electroproduction. The predictions described data over a large kinematic range in scale and in energy.
The dipole scattering amplitudes computed in this work are publicly available on the website https://hep.fjfi.cvut.cz/ along with instructions on how to use them.

VIII. ACKNOWLEDGEMENTS
We would like to thank Dionysios Triantafyllopoulos for fruitful discussions regarding