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

TL;DR

This paper presents a detailed solution to the impact-parameter dependent Balitsky-Kovchegov equation using a collinearly improved kernel, avoiding Coulomb tails and providing predictions consistent with experimental data.

Contribution

It introduces a novel solution to the impact-parameter dependent BK equation with an improved kernel that suppresses Coulomb tails and aligns well with experimental observations.

Findings

No Coulomb tails at large impact parameters.

Kernel suppresses large dipoles, affecting solution behavior.

Predictions agree with HERA and LHC data.

Abstract

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.

Equations

$$\frac{\partial N(\vec{r},\vec{b},Y)}{\partial Y}=\int d\overrightarrow{r_1}K(r,r_1,r_2)$$

$$\frac{\partial N(r,b,Y)}{\partial Y}=\int d\overrightarrow{r_1}K(r,r_1,r_2)$$

$$K_{\mathrm{LO}}(r,r_1,r_2)=\frac{{\alpha }_s^{\mathrm{nr}}}{2\pi }\frac{r^2}{r_1^2{r}_2^2},$$

$$K_{\mathrm{rc}}(r,r_1,r_2)=\frac{N_c{\alpha }_s(r^2)}{2{\pi }^2}\left(\frac{r^2}{r_1^2{r}_2^2}+\frac{1}{r_1^2}\left(\frac{{\alpha }_s(r_1^2)}{{\alpha }_s(r_2^2)}-1\right)+\frac{1}{r_2^2}\left(\frac{{\alpha }_s(r_2^2)}{{\alpha }_s(r_1^2)}-1\right)\right),$$

$$K_{\mathrm{rc}}^{\mathrm{bdep}}(r,r_1,r_2)=K_{\mathrm{rc}}(r,r_1,r_2)\Theta \left(\frac{1}{m^2}-r_1^2\right)\Theta \left(\frac{1}{m^2}-r_2^2\right),$$

$$K_{\mathrm{ci}}(r,r_1,r_2)=\frac{{\overline{\alpha }}_s}{2\pi }\frac{r^2}{r_1^2{r}_2^2}{\left[\frac{r^2}{\min (r_1^2,r_2^2)}\right]}^{\pm {\overline{\alpha }}_s{A}_1}{K}_{\mathrm{DLA}}(\sqrt{L_{r_{1}r}{L}_{r_{2}r}}),$$

$$K_{\mathrm{DLA}}(\rho )=\frac{J_1(2\sqrt{{\overline{\alpha }}_s{\rho }^2})}{\sqrt{{\overline{\alpha }}_s\rho }},$$

$$L_{r_{i}r}=\ln \left(\frac{r_i^2}{r^2}\right).$$

$${\overline{\alpha }}_s={\alpha }_s\frac{N_c}{\pi },$$

$${\alpha }_s={\alpha }_s(r_{\min }),$$

$${\alpha }_{s,n_f}(r^2)=\frac{4\pi }{{\beta }_{0,n_f}\ln \left(\frac{4C^2}{r^2{\Lambda }_{n_f}^2}\right)},$$

$${\beta }_{0,n_f}=11-\frac{2}{3}{n}_f.$$

$${\Lambda }_{n_f-1}=(m_f{)}^{1-\frac{{\beta }_{0,n_f}}{{\beta }_{0,n_f-1}}}({\Lambda }_{n_f}{)}^{\frac{{\beta }_{0,n_f}}{{\beta }_{0,n_f-1}}}.$$

$$r^2<\frac{4C^2}{m_f^2}.$$

$$N(r,b,Y=0)=1-\exp \left(-\frac{1}{2}\frac{Q_s^2}{4}{r}^{2}T(b_{q_1},b_{q_2})\right),$$

$$T(b_{q_1},b_{q_2})=\left[\exp \left(-\frac{b_{q_1}^2}{2B_G}\right)+\exp \left(-\frac{b_{q_2}^2}{2B_G}\right)\right].$$

$$K_{\mathrm{ci}}^1=\frac{{\overline{\alpha }}_s}{2\pi }\frac{r^2}{r_1^2{r}_2^2},$$

$$K_{\mathrm{ci}}^2={\left[\frac{r^2}{\min (r_1^2,r_2^2)}\right]}^{\pm {\overline{\alpha }}_s{A}_1},$$

$$K_{\mathrm{ci}}^3=K_{\mathrm{DLA}}(\sqrt{L_{r_{1}r}{L}_{r_{2}r}}).$$

$$K_{\mathrm{ci}}=K_{\mathrm{ci}}^1{K}_{\mathrm{ci}}^2{K}_{\mathrm{ci}}^3.$$

$$K_{\mathrm{rc}}^1=\frac{N_c{\alpha }_s(r^2)}{2{\pi }^2}\frac{r^2}{r_1^2{r}_2^2},$$

$$K_{\mathrm{rc}}^2=\frac{N_c{\alpha }_s(r^2)}{2{\pi }^2}\frac{1}{r_1^2}\left(\frac{{\alpha }_s(r_1^2)}{{\alpha }_s(r_2^2)}-1\right),$$

$$K_{\mathrm{rc}}^3=\frac{N_c{\alpha }_s(r^2)}{2{\pi }^2}\frac{1}{r_2^2}\left(\frac{{\alpha }_s(r_2^2)}{{\alpha }_s(r_1^2)}-1\right),$$

$$K_{\mathrm{rc}}=K_{\mathrm{rc}}^1+K_{\mathrm{rc}}^2+K_{\mathrm{rc}}^3.$$

$$N(r,Y)=\int d^2\vec{b}N(r,b,Y),$$

$$\mid\Psi_{T}^{i}(z,\vec{r},Q^{2})\mid^{2}=\frac{3\alpha_{\rm em}}{2\pi^{2}}e_{q_{i}}^{2}\Big{(}(z^{2}+(1-z)^{2})\epsilon^{2}K^{2}_{1}(\epsilon r)+m_{q_{i}}^{2}K^{2}_{0}(\epsilon r)\Big{)},$$

$$\mid\Psi_{L}^{i}(z,\vec{r},Q^{2})\mid^{2}=\frac{3\alpha_{\rm em}}{2\pi^{2}}e_{q_{i}}^{2}\Big{(}4Q^{2}z^{2}(1-z)^{2}K^{2}_{0}(\epsilon r)\Big{)}$$

$${ϵ}^2=z(1-z)Q^2+m_{q_i}^2,$$

$$\mid {\Psi }_{T,L}^i(z,\vec{r}){\mid }^2=\mid {\Psi }_T^i(z,\vec{r}){\mid }^2+\mid {\Psi }_L^i(z,\vec{r}){\mid }^2.$$

$$\frac{\mathrm{d}{\sigma }^{q\bar{q}}(\vec{r},x)}{\mathrm{d}\vec{b}}=2N(\vec{r},\vec{b},x).$$

$$\tilde{x}=x\left(1+\frac{4m_{q_i}^2}{Q^2}\right).$$

$$F_2(x,Q^2)=\frac{Q^2}{4{\pi }^2{\alpha }_{\mathrm{em}}}\int \sum _{i}d\vec{r}d\vec{b}dz\mid {\Psi }_{T,L}^i(z,\vec{r}){\mid }^2\frac{\mathrm{d}{\sigma }^{q\bar{q}}(\vec{r},\tilde{x})}{\mathrm{d}\vec{b}},$$

Full text

I Introduction

Evolution equations are powerful tools to study the high-energy, equivalently, small-$x$ limit of quantum chromodynamics (QCD) Altarelli (1982); Lipatov (1997); Gribov et al. (1983). The availability of quality data from HERA Newman and Wing (2014) and the LHC Akiba et al. (2016) as well as the need for reliable phenomenology for the proposal of new electron-ion facilities Accardi et al. (2016); Abelleira Fernandez et al. (2012) have given an extra impulse to the development of these tools.

In this work, the emphasis is placed on the Balitsky-Kovchegov (BK) evolution equation derived independently in the operator-product-expansion formalism by Balitsky Balitsky (1996), and by Kovchegov Kovchegov (1999, 2000) within the colour dipole approach Mueller (1990); Nikolaev and Zakharov (1991); Mueller (1994). It corresponds to the large-number-of-colours limit of the Jalilian-Marian-Iancu-McLerran-Weigert-Leonidov-Kovner (JIMWLK) evolution equations Jalilian-Marian et al. (1998a, b); Weigert (2002); Iancu et al. (2001a, b); Ferreiro et al. (2002). The BK equation describes the evolution with rapidity, $Y$, of the dipole-target scattering amplitude, $N(\vec{r},\vec{b},Y)$, where $\vec{r}$ is the transverse size of the dipole and $\vec{b}$ the impact parameter of the interaction.

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 Balitsky (2007); Kovchegov and Weigert (2007); Albacete and Kovchegov (2007). 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 Albacete et al. (2005, 2009, 2011); Cepila and Contreras (2015). A few other applications of these solutions are, for example, gluon production in heavy-ion collisions ALbacete and Dumitru (2010), single particle Lappi and Mäntysaari (2013) and $\mathrm{J/}\psi$ production in pp and pA collisions Ducloué et al. (2015) , di-hadron correlations in p–Pb interactions Albacete et al. (2019) and even the flux of atmospheric neutrinos Albacete et al. (2015); Bhattacharya et al. (2016).

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\equiv|\vec{b}|$, the so-called Coulomb tails, which generate an unphysical growth of the cross section Golec-Biernat and Stasto (2003). 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 Berger and Stasto (2011a). 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 Berger and Stasto (2011b). (A similar conclusion also holds for the solutions of the impact-parameter dependent JIMWLK equation Mäntysaari and Schenke (2018).) Nonetheless, this approach did a good job when confronted with HERA data on exclusive vector meson production Berger and Stasto (2013).

Recently, the kernel of the leading order equation has been improved by including the resummation of all double collinear logarithms Iancu et al. (2015a) as well as two classes of single logarithmic corrections Iancu et al. (2015b). 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 Cepila et al. (2019), 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 $\mathrm{J/}\psi$ at one energy without the need of any additional ad-hoc parameter or correction.

In this contribution the studies reported in Cepila et al. (2019) 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 $\phi$, $\mathrm{J/}\psi$, $\psi(2\mathrm{S})$, and $\Upsilon(\mathrm{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.

II Review of the formalism

II.1 The Balitsky-Kovchegov equation

The BK evolution equation reads Balitsky (2007); Kovchegov and Weigert (2007)

[TABLE]

where $r\equiv|\vec{r}\,|$, $r_{1}\equiv|\vec{r_{1}}|$, and $r_{2}\equiv|\vec{r_{2}}|\equiv|\vec{r}-\vec{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 are $b\equiv\vec{b}$, $b_{1}\equiv\vec{b_{1}}$, $b_{2}\equiv\vec{b_{2}}$. The kernel $K(r,r_{1},r_{2})$ is discussed below.

In this work, the solution to the BK equation is obtained under the assumption that the scattering amplitude $N(\vec{r},\vec{b},Y)$ depends solely on the sizes of the dipoles and of the impact parameter vectors. In practice, this means to solve the equation

[TABLE]

subjected to the condition that the angle between $\vec{r}$ and $\vec{b}$ is fixed. We chose to fix this angle at zero, meaning that these vectors are parallel.

II.2 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

[TABLE]

where the non-running coupling, $\alpha^{\rm nr}_{s}$, is fixed to a constant value.

The running coupling kernel $K_{{\rm rc}}(r,r_{1},r_{2})$ reads Balitsky (2007)

[TABLE]

where $N_{c}$ is the number of colours and $\alpha_{s}$ is the running coupling, which is further discussed in Sec. II.3.

The running coupling kernel with a cutoff to tame the Coulomb tails generated by the evolution in the impact parameter is given by Berger and Stasto (2011b)

[TABLE]

where $\Theta$ is the Heaviside function and $m$ a parameter to limit the size of daughter dipoles.

Finally, the collinearly improved kernel is Iancu et al. (2015b)

[TABLE]

where

[TABLE]

$J_{1}$ is the Bessel function (the inclusion of the Bessel function into the BK kernel has been previously discussed in Motyka and Stasto (2009)), the anomalous dimension is $A_{1}=11/12$, and

[TABLE]

The sign factor in the exponent $\pm\overline{\alpha}_{s}A_{1}$ takes the value of the plus sign when $r^{2}<\min(r_{1}^{2},r_{2}^{2})$ and the negative sign otherwise. For the running coupling

[TABLE]

the smallest dipole prescription is used throughout the computation according to

[TABLE]

where $r_{\min}=\min(r_{1},r_{2},r)$. This prescription was compared to other prescriptions in Iancu et al. (2015b), 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 Balitsky and Chirilli (2008).

II.3 Treatment of the coupling constant

In this work the running coupling is computed in the variable-number-of-flavours scheme, implemented according to

[TABLE]

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 Albacete et al. (2009). The variable $\beta_{0,n_{f}}$ is the leading order coefficient of the QCD beta-series and is given by relation

[TABLE]

The value of the QCD scale parameter $\Lambda^{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 relation Albacete et al. (2011)

[TABLE]

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 Z0 boson. In this way, $\Lambda_{5}$ is set with the use of the experimentally measured value of $\alpha_{s}(M_{Z})=0.1196\pm 0.0017$, where the Z0 mass is $M_{Z}=91.18$ GeV$/c^{2}$ Olive et al. (2014). 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

[TABLE]

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  Albacete et al. (2009). In this work, the coupling is frozen at $\alpha_{s}^{\rm sat}$ = 1 as in Iancu et al. (2015a).

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 $\alpha_{s}$ for two values of $C$: the one used here, $C=9$, and the one used in Iancu et al. (2015a), $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.

III Impact-parameter solution to the Balitsky-Kovchegov equation

III.1 Initial condition

The initial condition, already introduced in Cepila et al. (2019), 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 Golec-Biernat and Wusthoff (1998), 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 Kowalski and Teaney (2003); Watt and Kowalski (2008); Marquet (2007); Mäntysaari and Schenke (2016); Cepila et al. (2017). The main new ingredient with respect to the initial condition used in the previous studies Berger and Stasto (2011a, b, 2013) is the explicit separation of the contribution from the individual quark and anti-quark forming the dipole. The initial condition is given by

[TABLE]

where $b_{q_{i}}$ are the impact parameters of the quark and anti-quark forming the dipole and

[TABLE]

As a first attempt, the angle between $\vec{r}$ and $\vec{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 GeV2, 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\big{(}x_{0}/x\big{)}, 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 $\mathrm{J/}\psi$ off protons at a photon–proton centre-of-mass energy $\langle W_{\gamma p}\rangle=100$ GeV, where, as customary $x=(M^{2}+Q^{2})/(W_{\gamma p}^{2}+Q^{2})$ is used; here, $M$ represents the mass of the vector meson.

III.2 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 Cepila and Contreras (2015); Matas et al. (2016) was extended to the case of the impact-parameter dependent BK equation Cepila et al. (2019) and the solution is evolved in rapidity with a step of $\Delta Y$ = 0.01.

Fixed grids are used for $r$ and $b$. They are logarithmic grids of base 10 with 25 evenly-spaced points per order of magnitude, spanning the range from 10*-7* to 104 GeV*-1* for both the $r$ and $b$ variables. The integration over $\vec{r_{1}}$ is performed in polar coordinates, where $r_{1}$ is evaluated in the same grid as $r$ and the polar angle, denoted by $\theta_{rr_{1}}$, is evaluated in a fixed grid with 21 points separated by a constant step. The numerical integrations are performed applying Simpson’s method.

Since the transverse dipole vectors are related as $\vec{r}=\vec{r_{1}}+\vec{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 log10 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 $\vec{b_{1}}=\vec{b}+\vec{r_{2}}/2$ and $\vec{b_{2}}=\vec{b}-\vec{r_{1}}/2$ assuming a fixed angle between $\vec{r}$ and $\vec{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

IV.1 Behaviour of the collinearly improved kernel

As was shown in Cepila et al. (2019), 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 $\theta_{rr_{1}}=\pi/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

[TABLE]

The first term, $K_{\rm ci}^{1}$, is present already at the leading order if one considers a fixed value of the running coupling, $K_{\rm ci}^{2}$ takes into account the contribution from the single collinear logarithms, and $K_{\rm ci}^{3}$ resumms double collinear logarithms to all orders. The entire collinearly improved kernel is then given by the multiplication of all these factors as

[TABLE]

For the running coupling BK kernel, the separation in three parts is as follows:

[TABLE]

whereas the running coupling kernel is then given by the addition of these constituent terms as

[TABLE]

The contribution of the three terms is shown in Fig. 4 at $r$ = 1 GeV*-1* and $\theta_{rr_{1}}=\pi/2$ for each of the two kernels. The fact that the three terms are added in $K_{\rm rc}$, but multiplied in $K_{\rm ci}$ explains numerically the suppression. Even though, the first term is essentially the same for both kernels, the additive character of $K_{\rm 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_{\rm rc}^{2}$ and $K_{\rm rc}^{3}$ terms become dominant, whereas in the collinearly improved kernel, the $K_{\rm ci}^{1}$ term suppresses the total value.

The physical reason of this suppression can be traced back to the fact that large daughter dipoles do not follow the time-ordering prescription (that is, they would live longer than the parent dipole) built in when setting up the resummation that leads to the collinearly improved kernel Beuf (2014); Iancu et al. (2015b).

IV.2 Contribution of the kernel terms to the evolution

The suppression for large sizes of the daughter dipole in the kernel is translated as a suppression of the amplitude at large $b$ in the evolution. In this region only large $r_{1,2}$ 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}$) $\sim$ $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_{\rm 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_{\rm ci}^{3}$ term with the Bessel functions. This term originates from resumming double collinear logarithms. Note that also the term $K_{\rm ci}^{2}$, resumming single collinear logarithms, suppresses the large $b$ region.

IV.3 Behaviour of the solution to BK equation

The evolution of the scattering amplitude as a function of $r$ for different fixed values of $b$ is shown in the upper panels of Fig. 6, while the lower panels of the same figure show the evolution as a function of $b$ for two fixed values of $r$. A two-dimensional view of the amplitude at two stages of the evolution is shown in Fig. 7. The amplitude decreases fast for small dipole sizes as expected. The suppression of large dipole sizes imposed in the initial condition is lifted with evolution. Similar behaviour was observed in previous studies, e.g. Berger and Stasto (2011a). 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 Contreras et al. (2019).

Finally, Fig. 8 shows $N(r,Y)$, defined as

[TABLE]

for different dipole sizes and for two kernels, the running coupling and the collinearly improved. 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 GeV2 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 Cepila and Contreras (2015).

Another interesting observation is that $N(r,Y)$ is related to the $\sigma_{0}$ parameter introduced in studies based on the rcBK equation without impact-parameter dependence. Basically, $\sigma_{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 Albacete et al. (2011). Figure 8 justifies the order of magnitude of these values from the perspective of an impact-parameter dependent computation.

V Deep Inelastic Scattering

V.1 Structure function and reduced cross section

Due to the fact, that the dipole lives much longer than the typical interaction time, the computation of the total deep-inelastic scattering (DIS) cross section can be written as the convolution of separate terms. One of them is the wave function representing the probability of a virtual photon splitting into a quark-antiquark dipole. Here formulas and notation of Golec-Biernat and Wusthoff (1998) are used:

[TABLE]

and

[TABLE]

for the transverse and longitudinal polarisation of the incoming photon, respectively, where $z$ is the fraction of the total longitudinal momentum of the photon carried by the quark, $K_{0}$ and $K_{1}$ are the MacDonald functions, $Q^{2}$ is the virtuality of the probing photon, $e_{qi}$ is the fractional charge (in units of elementary charge) of quark $i$, $\alpha_{\rm em}$ = 1/137 and

[TABLE]

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 Iancu et al. (2015b)); 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

[TABLE]

According to the optical theorem, one can link the dipole–target cross section to the scattering amplitude by

[TABLE]

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 Golec-Biernat and Wusthoff (1998)

[TABLE]

Using these ingredients, the relation for the computation of the structure function in the dipole model framework is

[TABLE]

and the reduced cross section is computed as

[TABLE]

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

[TABLE]

V.2 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 Aaron et al. (2010a) 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

[TABLE]

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 Aaron et al. (2010a).

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 Balitsky and Chirilli (2008), and a large effort is being done to use it for phenomenology Lappi and Mäntysaari (2016); Beuf (2014, 2017); Ducloué et al. (2017); Hänninen et al. (2018). There are also recent developments regarding the most adequate variable to evolve the scattering amplitude Ducloué et al. (2019).

VI Production of vector mesons

VI.1 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 $\mathrm{V}$ is given by

[TABLE]

where $\mathcal{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., Kowalski et al. (2006)) and takes the following form

[TABLE]

where the subscripts $T$, $L$ denote the contribution from the virtual photon with transverse, respectively longitudinal, polarisation, $\Psi_{\gamma^{*}}$ is the wave function of a virtual photon which fluctuates into a dipole, $\Psi_{\mathrm{V}}$ represents the wave function of the vector meson, and $\vec{\Delta}^{2}\equiv-t$, the square of the four momentum transferred in the proton vertex. Under the assumption of large photon-proton centre-of-mass energies $W_{\gamma\mathrm{p}}$,

[TABLE]

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 $q\bar{q}$ 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

[TABLE]

and

[TABLE]

with $\hat{e}_{f}$ being the effective charge of the given vector meson, $\epsilon$ defined by Eq. (28), and the parameter $\delta$ is a switch to include ($\delta=1$) or exclude ($\delta=0$) the non-local term in the longitudinal contribution. The scalar part $\phi_{T,L}$ of the wave function is, in general, model dependent. For our studies, we use the boosted Gaussian model Nemchik et al. (1994, 1997); Forshaw et al. (2004) in which the $\delta$ 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 Bendova et al. (2019).

The total exclusive cross section to produce a vector meson is given by the sum of the transverse and the longitudinal contributions defined by Eq. (36). Moreover two important corrections have to be applied. The derivation of the formula for the exclusive vector meson cross section is performed under the assumption that the scattering amplitude $\mathcal{A}_{T,L}(x,Q^{2},\vec{\Delta})$ is purely imaginary. The real part of the amplitude can be accounted for by the extra term $(1+\beta^{2})$ in Eq. (37), where $\beta$ is the ratio of real to imaginary parts of the scattering amplitude, for details see Kowalski et al. (2006). 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 off-diagonal gluon distribution for vector meson production. This effect can be accounted for by multiplying the scattering amplitude by a factor $R_{g}^{T,L}$, called the skewedness correction Shuvaev et al. (1999).

VI.2 Comparison to data

Using the model described in this paper, the cross sections for exclusive photo- and electroproduction of $\phi$, $\mathrm{J/}\psi$, $\psi(2\mathrm{S})$, and $\Upsilon(\mathrm{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 a comparison of our predictions for the $|t|$-distributions and the total cross sections with HERA H1 Aaron et al. (2010b) and ZEUS Chekanov et al. (2005) data for the exclusive electroproduction of the $\phi$ meson for several values of $Q^{2}$ is shown. The predictions give a very good description of the available data, especially at low photon virtualities.

The predictions for the exclusive production of the $\mathrm{J/}\psi$ meson are compared with the experimental data from H1 Aktas et al. (2006); Alexa et al. (2013) and ALICE Abelev et al. (2014); Acharya et al. (2019) experiments in Figs. 14-16, for several different measurements of kinematic observables. In the left panel of the Fig. 14, the comparison of the $|t|$-distribution of the photoproduction cross section is presented. The predictions give very good agreement with the data at energies $W_{\gamma\mathrm{p}}=55\;\mathrm{GeV}$ and $W_{\gamma\mathrm{p}}=100\;\mathrm{GeV}$. The result for $W_{\gamma\mathrm{p}}=78\;\mathrm{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_{\gamma\mathrm{p}}=100\;\mathrm{GeV}$. Since the value of $W_{\gamma\mathrm{p}}$ from the experimental data is a mean value estimated from a measured energy range, the result of the model can be considered satisfactory. 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_{\gamma\mathrm{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_{\gamma\mathrm{p}}$, however at larger values ($\sim 10^{2}\;\mathrm{GeV}$), the predictions are underestimated when compared to experimental photoproduction data. We have also obtained total cross section for the $\mathrm{J/}\psi$ 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_{\gamma\mathrm{p}}$, however at high energies the result is again underestimated when compared to data.

Also, the exclusive cross section of the $\psi(2\mathrm{S})$ meson was calculated within the model. The experimental data are not available for the total cross sections, but only for a ratio of the $\psi(2\mathrm{S})$ to $\mathrm{J/}\psi$ cross sections, the predictions for these ratios for photoproduction and electroproduction at $Q^{2}=16\;\mathrm{GeV^{2}}$ are calculated and compared to data from H1 Adloff et al. (2002), and ZEUS Abramowicz et al. (2016), 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 $\Upsilon(\mathrm{1S})$ meson is presented in Fig. 17. The prediction is compared with experimental data obtained at HERA by H1 Adloff et al. (2000) and ZEUS Chekanov et al. (2009) experiments. It is also compared with the two latest measurements – in proton-proton collisions at $\sqrt{s}=7\;\mathrm{TeV}$ and $\sqrt{s}=8\;\mathrm{TeV}$ by LHCb Aaij et al. (2015), and in proton-lead collisions at $\sqrt{s}=5.02\;\mathrm{TeV}$ by the CMS experiment Sirunyan et al. (2019). 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 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 this paper. Our work has been partially supported from grant LTC17038 of the INTER-EXCELLENCE program at the Ministry of Education, Youth and Sports of the Czech Republic, by grant 17-04505S of the Czech Science Foundation, GA ČR and the COST Action CA15213 THOR. Computational resources were provided by the CESNET LM2015042 grant and the CERIT Scientific Cloud LM2015085, provided under the program Projects of Large Research, Development, and Innovations Infrastructures.