Markov Chain Mote Carlo solution of BK equation through Newton-Kantorovich method

We propose a new method for Monte Carlo solution of non-linear integral equations by combining the Newton-Kantorovich method for solving non-linear equations with the Markov Chain Monte Carlo (MCMC) method for solving linear equations. The Newton-Kantorovich method allows to express the non-linear equation as a system of the linear equations which then can be treated by the MCMC (random walk) algorithm. We apply this method to the Balitsky-Kovchegov (BK) equation describing evolution of gluon density at low x. Results of numerical computations show that the MCMC method is both precise and efficient. The presented algorithm may be particularly suited for solving more complicated and higher-dimensional non-linear integral equation, for which traditional methods become unfeasible.


Introduction
The Large Hadron Collider (LHC) at CERN provides an opportunity to scan parton densities in the proton over a wide domain of parton kinematics. This allows for detailed studies of dynamical effects taking place during evolutions of partons. An example of the dynamical phenomena which is particularly interesting in hadronic processes is saturation of gluon density [1]. At high energies the dominant contribution to evolution of system of partons comes from splittings of gluons and this leads to rapid growth of the gluon density and, as a consequence, to fast rise of corresponding cross sections. The unitarity constraints suggest that eventually the growth of the gluon density should slow down due to possible effects of fusion of the gluons, leading to its saturation. And indeed, there is a growing evidence that the saturation occurs in high-energy hadron collision processes [2][3][4][5].
The physics of the saturation at an inclusive level is described within the perturbative QCD by [6][7][8][9][10][11], and at an exclusive level by the equations proposed in [12][13][14]. The standard approach in search of the saturation with the JIMWLK/BK evolution equation is to solve the equation that provides the gluon density, and then to convolute the solutions with appropriate matrix elements which specify the final states. This approach has some limitations, since it does not allow for the full simulation of a scattering event, as is for example modelled by Monte Carlo event generators [15][16][17][18][19]. The Monte Carlo event generators allow for exact treatment of kinematical effects, storing information on emitted partons, etc. The particularly useful method of performing Monte Carlo simulation is based on a Markov Chain (random walk) approach [20]. In this approach, the evolution process occurs over evolution 'time' which could be, for example, an energy scale. Such an evolution can be interpreted as a Markovian probabilistic process. The main advantage of this approach is that one performs the full Monte Carlo simulation in a forward process, without the need to pretabulate the solution of the considered equation. In the so-called backward evolution method, first the appropriate equation is solved and the corresponding parton density is pretabulated, and then the actual Monte Carlo evolution (random walk) is performed backward in a 'time' variable to simulate the scattering process with a probability distribution JHEP07(2013)097 given by the respective parton density. The application of the Markov Chain Monte Carlo (MCMC) algorithm is, however, not straightforward for equations which model the saturation effects, since it works only for the linear evolution equations. To our best knowledge, such an algorithm has not been, so far, applied to the non-linear evolution equations.
In the present paper we develop a method which allows to perform the MCMC-based solution of the non-linear equation. 1 We apply this method to the BK equation. The key idea is to apply a well-convergent method for solving the non-linear integral equation and to combine it with a Monte Carlo algorithm designed for solving the linear integral equations. We have found that particularly well-suited for this purpose is the Newton-Kantorovich method [22]. It relies on representing the non-linear integral equation as a system of the linear equations (see eq. (2.7)-(2.9)), to which the MCMC algorithm can be applied. The whole procedure can be done in iterative manner and it does not require one to provide the solution of the considered equation in advance. Furthermore, it can be used for solving the exclusive saturation equations, as proposed in [12][13][14], and even more complicated and higher-dimensional non-linear integral equations, for which other numerical methods are unfeasible (inefficient, ustable, etc.). This might also be a first step in constructing the Monte Carlo event generator for modelling the saturation effects in the fully exclusive way.
The paper is organized as follows. In section 2 we introduce the Newton-Kantorovich method for the BK equation. In section 3 we describe the Markov Chain Monte Carlo algorithm. In section 4 we combine the Newton-Kantorovich method with the Monte Carlo algorithm to provide the solution of the BK equation. We also compare our solution with the one provided by the BKSolver package [23].
where Φ 0 (x, k 2 ) is a driving term,ᾱ s = (N c α s )/π (in our calculations we use α s = 0.2), k ≡ k ⊥ is the transverse gluon momentum, x is the fraction of the longitudinal proton momentum carried by the gluon, and hereinafter we set R = 1/ √ π. First, we perform the following change of variables: y = − ln x, t = y + ln z ⇒ x = e −y , x/z = e −t and simplify the notation by skipping exponents of arguments of the above functions, i.e. Φ(e −y , . . .) → Φ(y, . . .), etc., to obtain:

JHEP07(2013)097
Introducing a dimensionful constant µ 2 , the dimensional integration variable l 2 can be replaced by λ = ln(l 2 /µ 2 ), for which we have dλ = dl 2 /l 2 . Introducing also κ = ln(k 2 /µ 2 ), we get where the notation is again simplified: we use κ instead of k 2 in arguments of Φ and Φ 0 , and drop the dependence on the scale µ 2 , which is obviously hidden in both functions. This is the two-dimesional non-linear integral equation of the form: One can linearize the above equation by expanding the kernel K (y, t, κ, λ, Φ(t, λ)) in the Taylor series with respect to Φ(t, λ) about someΦ(t, λ) and retaining only the first two terms: Introducing the function Ψ(t, λ) = Φ(t, λ) −Φ(t, λ) and assuming |Ψ(t, λ)| ≪ 1, we can replace eq. (2.3) with the following set of equations: (2.10)

JHEP07(2013)097
The above set of equations can be solved by iteration, which leads to the Newton-Kantorovich form of the BK equation: Φ n (y, κ) = Φ n−1 (y, κ) + Ψ n−1 (y, κ), (2.11) As one can see, instead of the single non-linear integral equation (2.4) we have now the iterative series of the linear integral equations (2.12), associated with the auxiliary integrals of eq. (2.13). This can be solved by the standard iteration (successive approximation) method. The two-dimesional integrations can be performed directly with the standard numerical quadratures or, alternatively, one may expand the integrands in series of the Chebyshev polynomials, at least in one integration variable. The main advantage of the above decomposition is that the integral equation (2.12) is linear, and thus one can try to solve it by using the MCMC algorithm.

MCMC method
Our goal in this section is to construct a MCMC solution of eq. (2.12) which is the Volterra-Fredholm linear integral equation of the second kind. We can write immediately its iterative solution: Ψ n−1 (y, κ) = Λ n−1 (y, κ) (3.1) Since the integration limits do not depend on the variable κ, there is no ordering in the integration variable λ and at any step it can take an arbitrary value. Due to the ordering in the integration variable t, it will play a role of the evolution time in the corresponding MCMC algorithm. We propose the following MCMC algorithm: 1. Start a random walk (Markov chain) from the point (t 0 , λ 0 ) = (y, κ).

JHEP07(2013)097
where η τ (ξ) is the pdf of the variable ξ depending on the parameter τ (if it does not depend on this parameter, then ξ can be generated completely independently of τ ).

To each trajectory
2) assign the von Neumann-Ulam weight 2 [20]: where is the probalility of a single jump beyond y 0 from the point t.
Instead of the von Neumann-Ulam weight one may use the Wasow weight 3 [24]: 5. Repeat the above steps N times and compute the MCMC estimate of Ψ n−1 (y, κ) as well as its statistical error (standard deviation): where w k (y, κ) is the trajectory weight (of eq. (3.3) or eq. (3.5)) computed in the kth repetition of the above steps 1-4. One can prove that expectation values of the weights w(y, κ) of eqs. (3.3) and (3.5) satisfy the equation (2.12). The most straightforward way to do this is to first obtain general expressions for contributions to the weights coming from the trajectory of the length m, and then, based on that, construct the corresponding expectation values.

JHEP07(2013)097 4 Numerical results
In this section we present an implementation of the MCMC algorithm. We perform computations on a 2-dimensional lattice of points -in the rapidity y and in the dimensionless variable κ, corresponding to the transverse momentum k ⊥ . The results presented here correspond to the lattice with 100 points in y distributed linearly from 0.0 to 8.1, and 128 points in κ linearly spread in the range [0.0, 10.6]. The k ⊥ dimension is introduced to the problem through the constant µ 2 which shows up in the driving term Φ 0 (y, κ): In our computations we have used µ 2 = 5 · 10 −3 GeV 2 . In eq. (2.3) one can see that the integration over λ goes to infinity. In order to perform numerical calculations we need to introduce a certain cut-off. The driving term of eq. (4.1) as well as the solution of the BK equation vanish for large k ⊥ , thus introducing the upper cut-off on κ does not affect the solution considerably.
For the pdfs ρ(τ ) and η(ξ) we use the exponential distributions: and thus the random variables τ i and λ i can be generated as follows: where U i and V i are the random variables uniformly distributed between 0 and 1, i.e. U i , V i ∈ U (0, 1). This choice does its job in the case of the above BK equation, however one can improve the convergence of the MCMC method by using the pdfs that are better adjusted to the problem. Ideally, the product of these pdfs should be as close as possible to the kernel K ′ Φ , so that all the weights v in eq. (3.3) be close to 1. In fact our choice of the pdfs seems to be good enough as we have reached a sufficient precision generating only 1000 trajectories for each iteration of eq. (2.12). The results presented here correspond to 15 iterations of the set of equations (2.11)-(2.13). Without special optimisations it took only about 20 minutes of CPU time to generate all the results on a 2.2 GHz Intel Pentium Dual-Core processor with the GNU/Linux operating system using only one CPU core.
As stated in the previous section, one can use either the von Neumann-Ulam or the Wasow weights in the MCMC procedure. Our implementation of the MCMC algorithm has been tested with both of them, giving the same results (differences not visible in the plots like the one presented below).
In figure 1 we show the results of the numerical solution of the BK equation. One can see the profile plots of our MCMC solutions together with the reference solutions obtained with BKsolver. The latter program evolves the solution of the BK equation in rapidity based on the differential version of the equation. We have run BKsolver with the same parameter ranges as stated above and required the same number of points in the output lattice. The plots presented here have been then obtained using two-dimensional bilinear interpolation. In figure 1(a) we present the solutions in the k ⊥ -profile for three different rapidity values: y = 2, 5, 8. For each of them the results from our MCMC algorithm and from BKsolver are shown. In the lower part of the plot one can see the relative difference between the MCMC solution and the reference BKsolver one. As one can see, these two solutions agree at the level below 0.1 %. Similarly, figure 1(b) contains three y-profiles, each for different k ⊥ values: k ⊥ = 0.1, 1, 10 GeV. The results from both the MCMC implementation and BKsolver are shown as well as their relative difference. The agreement between the two solutions is again below 0.1 %.
One might have realized that the results in figure 1 are shown in the narrower rapidity and k ⊥ ranges than given at the beginning of this section. We have simply skipped some points on the lattice boundaries where the agreement between the two methods is slightly worse. It is caused by such factors as the interpolation errors and/or the finite number of the lattice points, rather than by a problem of the MCMC algorithm itself. Small fluctuations of the relative differences in both plots (a) and (b) are due to finite numbers of points in the lattices and approximations of the interpolation procedure.
Generally, with these results we have proved that the MCMC algorithm is applicable for solving the BK equation and, indeed, it gives good numerical results.

Summary and outlook
In this paper we have developed a general method to solve the two-dimensional non-linear integral equation via Monte Carlo techniques. Our method relies on combining the robust Newton-Kantorovich procedure for solving the non-linear integral equations with the Markov Chain Monte Carlo algorithm. The method is powerful and can be applied to solving complicated, high-dimensional non-linear integral equations, where the traditional JHEP07(2013)097 methods become inefficient. It can also open a window to construction of a Monte Carlo event generator based on the non-linear integral equations, which will allow to study saturation effects in the fully exclusive processes.
We have applied the MCMC algorithm to the BK equation and compared the results with the ones obtained by using the traditional methods, i.e. the solution of the BK equation as an integro-differential equation, implemented in the BKsolver package. The agreement within 0.1% have been found. The presented MCMC algorithm is general, and thus it can also be applied to the exclusive form of the BK [12] and KGBJS [12] evolution equations. This we leave, however, for the future studies.