Some numerical studies of the evolution of generalized parton distributions

We study the evolution behavior of generalized parton distributions at small longitudinal momentum fraction. Particular attention is paid to the ratio of a generalized parton distribution and its forward limit, to the mixing between quarks and gluons, and to the dependence on the squared momentum transfer t.


Introduction
A characteristic property of generalized parton distributions (GPDs) is their renormalization scale dependence, described by evolution equations whose derivation led to the very discovery of these functions more than a decade ago [1]. On a practical level, the scale dependence of GPDs is of direct importance for the quantitative description of exclusive scattering processes. Moreover, understanding general features of the evolution behavior should be helpful for developing realistic models and parameterizations of GPDs. The question how a given input distribution changes when evolved to higher scales has been addressed in several studies, both numerically and analytically [2][3][4][5][6][7][8]. Further progress has been achieved recently [9,10] by constructing explicit solutions of the evolution equations with methods that generalize the familiar Mellin moment inversion for parton density functions (PDFs).
The aim of the present contribution is to study a number of aspects in the evolution of GPDs at a numerical level. We will largely concentrate on the value of the GPDs at x = ξ, which at leading order in α s determines the imaginary part of scattering amplitudes, and via dispersion relations also gives their real part up to a ξ independent constant [11]. Furthermore we will focus on the region of small ξ, where the behavior of distributions can be conveniently approximated by a power-law behavior. We will pay special attention to the mixing between the gluon GPD H g (x, ξ, t) and H S (x, ξ, t) = n f q H q (x, ξ, t) − H q (−x, ξ, t) , (1) whose forward limit is the familiar singlet combination of quark and antiquark PDFs. For comparison we will also consider H u−d (x, ξ, t) = H u (x, ξ, t) − H d (x, ξ, t) as a representative of the non-singlet sector. After specifying in Sect. 2 the GPD model used as initial condition for the evolution, we devote most of Sect. 3 to a quantitative study of the old question how the ratio of GPDs and PDFs behaves when evolved to higher scales. We shall in addition take a look at the behavior of the GPDs around x = ξ. In Sects. 4 and 5 we turn to the dependence of GPDs on the squared momentum transfer t. Both theoretical considerations [12] and lattice QCD calculations [13] indicate that this dependence is correlated with the one on the longitudinal variables x and ξ. Since evolution affects the x dependence at given ξ and t, it also affects the t dependence at given x and ξ in a nontrivial fashion, which we will quantify in two model scenarios.
For our calculations we have used the numerical code of [14], which provides a numerically fast and stable implementation of GPD evolution at leading order (LO) in α s . The effects of next-to-leading (NLO) and next-to-next-to-leading order (NNLO) terms in the evolution kernels have been studied [8,15] and are known to be important, especially at small ξ in the gluon and singlet sector. This should be kept in mind as a caveat when interpreting our results, but we think that a study at LO is still of some relevance. On one hand, the arguments in [5,7] about the pattern of evolution to higher scales are based on the LO kernels, so that this order is adequate to test the numerical validity of these arguments. On the other hand, evolution effects on the t dependence are barely known at all, and LO results should at least provide a valid starting point for further investigation.
In the evolution kernels and the running coupling we take n f = 4 for m c ≤ µ < m b and n f = 5 for µ ≥ m b , with the charm and bottom quark masses m c = 1.3 GeV and m b = 4.5 GeV used in the CTEQ6 parton analysis [16], which we use for calculating the model GPDs at the starting scale of evolution. We furthermore follow the CTEQ6 analysis in taking the twoloop running coupling with Λ (4) = 326 MeV and Λ (5) = 226 MeV, which corresponds to α (4) s (m c ) = 0.40 and α (5) s (M Z ) = 0.118. We shall not consider scales below µ = 1.3 GeV, which we regard as a compromise between starting evolution at a "low scale" and staying in a region where α s is not so large that the LO approximation becomes more and more questionable.

Initial conditions
At the starting scale of evolution, we use the Musatov-Radyushkin ansatz, which is based on double distributions [7]. With the conventional definitions of H q and H g , given e.g. in [17], we can write this ansatz as In this work we will use different values of the profile parameter b, which for simplicity will always be taken equal for all quark and gluon distributions. The ansatz (3) has been extensively used in the literature so far. One should keep in mind that it does not exhaust the possibilities of modeling, and other approaches [18][19][20][21] are being pursued in the literature. As we will see, this model does however provide enough flexibility to address a number of important questions. The model also permits useful analytic approximations at small ξ. At x = ξ the integrals in (3) are restricted to β < 2ξ, so that for ξ ≪ 1 one can neglect the β dependence in h b (β, α). Approximating 1 + ξ by 1 in the integration limits, one then has [7] H i (ξ, ξ, t) with i = q, g. We will use this approximation shortly.

Evolution at fixed t
In this section we study the evolution of GPDs at a fixed value of t. We take t = 0 and do not display this variable for brevity. To quantify the difference between generalized and usual parton distributions we use the conventional skewness ratios where we have explicitly displayed the dependence on the scale µ in the distributions.
As is well known, the PDFs obtained from fits to data follow an approximate power-law behavior at small x, at given µ, where a and λ depend of course on the parton species. With the ansatz (3) for GPDs this leads to a power-law behavior of the GPDs at small ξ according to (5), with the same powers λ as for the corresponding PDFs. The skewness ratios at small ξ are readily obtained as at the scale µ where the ansatz (3) is made. Numerically, we find that the approximate power-laws (7) and (8) remain valid under evolution to higher scales, with powers λ that depend on µ but remain the same for the forward distributions and the GPDs. Based on the considerations using the Shuvaev transformation, it has been argued in [5] that at small ξ and high enough scale, the skewness ratio should be given by [4,22] R g for gluons and quarks, respectively. Here λ is the power in (7) at the scale where R g (ξ, µ) or R q (ξ, µ) is evaluated. More precisely, the ratios in (10) are obtained if (7) holds and if all Gegenbauer moments of the GPD in question are independent of ξ. Musatov and Radyushkin [7] have shown that at small x and ξ this condition is tantamount to the GPD being given by (3) with b = λ + 1, for both gluon and quark distributions. Indeed, one can easily check that for i = g, q. Using a different line of arguments, the authors of [10] also expect that (10) should become valid after LO evolution to high scales, provided that one takes a particular joint limit of large µ and 1/ξ. The relations (10) are often used to calculate highenergy scattering amplitudes, so that it is important to test under which conditions they may be assumed to hold. We have taken the double distribution model (3) with the CTEQ6L distributions [16] at µ 0 = 1.3 GeV as input. After LO evolution to a scale µ, we have fitted effective power laws for g(x) and H g (ξ, ξ), and we have evaluated the skewness ratio R g (ξ, µ) from (6). In analogy we have determined power laws and ratios R S and R u−d for the combinations H S and H u−d introduced in Sect. 1.
Let us first discuss the power-law behavior (7) of the PDFs, which is not exact and only valid in a certain range of x. We fitted power-laws to the CTEQ6L parameterization for g(x), S(x) and u(x) − d(x) in the three intervals [10 −5 , 10 −4 ], [10 −4 , 10 −3 ] and [10 −3 , 10 −2 ]. The resulting powers for the gluon and quark singlet distributions are shown in Fig. 1. We see a clear x dependence of the effective power λ, especially at larger x. In Fig. 2 we show the powers obtained in the interval 10 −4 < x < 10 −3 for a larger range of µ. We note that under evolution the powers for the gluon and the quark singlet become similar but remain different up to very high µ. This effect has already been pointed out in [23]. For the non-singlet distribution u − d the effective power λ is between −0.41 and −0.42 in all three x intervals. It changes by less than 1% under evolution in the µ range corresponding to Figs. 1 and 2.
According to (5) there is no simple relation between the ranges of x and ξ in which the same power-law behavior should approximately hold for a PDF and the corresponding GPD. For simplicity we have fitted H g (ξ, ξ), H S (ξ, ξ) and H u−d (ξ, ξ) to power laws (8) in the same ξ intervals that we took for the PDFs. An example of such a fit is shown in Fig. 3, where we see that H g (ξ, ξ) indeed follows an approximate power-law over about one order of magnitude in ξ but not over the full range of the plot. We find that corresponding powers λ for PDFs and GPDs differ by at most 3% in the respective x and µ ranges of Figs. 1 and 2. An exception is the quark singlet distribution in the interval 10 −3 < x < 10 −2 , where the power for the GPD is higher than that for the PDF by 5% to 10%. This is not surprising, given that already in Fig. 1 we see a more rapid change of the effective power at higher x. For definiteness we will evaluate R b (λ) and R Sh (λ) with the powers fitted to the PDFs. We have checked that our conclusions do not change when taking the powers for the GPDs instead.
We note that in a specific joint limit of large µ and 1/ξ, the solutions of the LO evolution equations for PDFs exhibit so-called double logarithmic scaling [24]. In this case one obtains The effective powers in (7) are then larger for the gluon than for the quark singlet distribution and depend logarithmically on x. Double logarithmic scaling for GPDs in the region x ≥ ξ has been discussed in [10].
We have evaluated the skewness ratios R(ξ, µ) from the evolved GPDs and PDFs for ξ = 3.2 × 10 −5 , 3.2×10 −4 and 3.2×10 −3 . This is compared with R b (λ) and R Sh (λ) calculated with λ from our fits of the PDFs at the corresponding scale µ and in the corresponding At the starting scale µ 0 = 1.3 GeV the curves for R g (ξ, µ) and R g 2 (λ) coincide as they should, whereas for increasing µ they become different. This means that one obtains different results for H g (ξ, ξ; µ) when making the ansatz (3) at scale µ or when making it at scale µ 0 and then evolving the GPD. The difference is however fairly small.
The curves for R g (ξ, µ) and those for R g Sh (λ) = R g λ+1 (λ) in Figs. 4 and 5 are rather close to each other. At the starting scale they hardly differ at all, which reflects a particularity of the model ansatz (3) for gluons. This is because the ratio R g b (λ) in (9) has a very weak b dependence for small λ: varying b from 1 to ∞ one obtains for instance R g b (0.1) between 1.072 and 1.088. With increasing λ the b dependence grows only slowly, with R g b (0.3) between 1.231 and 1.307 for b between 1 and ∞. This is also seen in the left panel of Fig. 5, where R g b (λ) is given for several b values. The solid curve in this figure is for b = 2 in the initial condition, but the corresponding results for b = 1 or b = 8 differ by less than 0.5% in the µ range of the figure.
Obviously it is hard to see whether R g (ξ, µ) tends to R g Sh (λ) under evolution if the two functions are already close at the starting scale. To investigate this further    we take a variant of (3), namely with h b (β, α) as in (4). This corresponds to a double distribution representation for , and one readily verifies that it gives Mellin moments of H g (x, ξ, t) with a polynomial dependence on ξ as required by Lorentz invariance. An analogous representation was first discussed for the quark GPD of the pion [25] and was recently found to be relevant for polarized gluon GPDs [26]. In the case of H g the ansatz (12) has the peculiar property of giving a zero at x = 0 that is not required by symmetry and quickly disappears under evolution. One may therefore not take this model too seriously, but it serves the purpose of giving a skewness ratio sufficiently different from the one obtained with the more conventional ansatz (3). This is shown in the right panel of Fig. 5, where the dot-dashed curve corresponds to initial conditions (12) for H g and (3) for H S , with b = 2 in both cases. We see that the ratio R g in the two models indeed tends to a common value after evolution. This value it not exactly equal to R g Sh (λ) but differs from it by less than 2%. Such a small difference should not be regarded as significant: the form (10) of R g Sh (λ) is obtained in [4,22] from an integral of g(x) over x from ξ/2 to 1, assuming the power behavior (7) in the entire interval. This is clearly an approximation.
We now turn to the skewness ratio for the quark singlet distribution, which is shown in Fig. 6. In contrast to the gluon case, different values of b in the ansatz (3) lead to significantly different skewness ratios at the starting scale. Evolution to higher µ brings the curves of R S (ξ, µ) for different initial conditions closer     to each other. As in the case of gluons, they do not exactly approach the curve we calculate for R S Sh (λ), but again this should not be regarded as significant since the power-law (7) with a fixed value of λ is only an approximation for a certain x range.
In the left panel of Fig. 6 we also see a clear difference between the evolved ratios R S (ξ, µ) and the curves for R S b (λ) with λ taken at the corresponding scale µ. In general there is hence a notable dependence of H S on the scale where the ansatz (3) is made, especially for larger values of b. We find the dependence less pronounced for b = 1, in agreement with what was found in [7] for ξ = 5.26 × 10 −2 .
Based on the inversion of Gegenbauer moments, it was argued in [10] that quark distributions H q (x, ξ) should develop a singular derivative (∂/∂x)H q (x, ξ) at x = ξ after evolution. To investigate this, we have numerically calculated (∂/∂x)H S (x, ξ) from the difference quotient for successive points in x, which around x = ξ were spaced in intervals of 9 × 10 −7 . In Fig. 7 we plot H S (x, ξ) together with its logarithmic derivative (∂/∂x) ln H S (x, ξ). Taking the derivative of (3) and making the same approximations which lead to (5), one finds that (∂/∂x)H S (x, ξ) is singular at x = ξ for b ≤ 1+λ. Indeed, we see in the figure that at the starting scale the derivative has a singularity for b = 1 but remains finite for b = 2. Under evolution a singularity develops for b = 2, whereas for b = 1 the logarithmic derivative hardly changes. Notice that in the curves for H S (x, ξ) one can barely recognize that the tangent at x = ξ should be vertical: this illustrates the limitations of rendering a weakly singular derivative in a plot. In contrast to the quark case, the ansatz (3) with b = 1 or b = 2 gives a finite value of (∂/∂x)H g (x, ξ) at x = ξ. We have checked numerically that for both initial conditions the derivative remains finite under evolution, in agreement with what one expects from the analytical representation in [10].
To conclude this section, we briefly investigate the quark non-singlet distribution H u−d . In Fig. 8 we show the skewness ratio for different values of b in the initial condition. The corresponding curves for R u−d b (λ) are not shown: they coincide with those for R u−d (ξ, µ) at the starting scale and then remain essentially flat since the effective power λ hardly changes with µ in this case. Under evolution to high scales the curves for different b approach each other and the one for R u−d Sh (λ), although much more slowly than for R g or R S .

Ansatz for the t dependence
To investigate the change of the t dependence with evolution, we will use the model (3) with b = 2 for quarks and gluons. We thus need an ansatz for the GPDs at zero skewness ξ but finite t, which is described in this section. In all cases we assume an exponential t dependence that is correlated with x. For the valencetype combination of GPDs we take the form proposed in Ref. [27]: The values α ′ v = 0.9 GeV −2 , B u = B d = 0.59 GeV −2 , A u = 1.22 GeV −2 and A d = 2.59 GeV −2 together with the CTEQ6M parameterization for q v (x) at µ = 2 GeV lead to a good description of the data for the Dirac form factors F 1 (t) of proton and neutron, which are obtained by combining dx H q v (x, 0, t) for u and d quarks with the appropriate charge factors.
For small x we can approximate (14) as f q (x) ≈ α ′ v ln(1/x) + B q and thus have where in the second step we have assumed a small-x behavior of the valence quark distributions as in (7).
Since the x dependence of (15) is a power-law, the integral in (5) can be performed as in Sect. 3, and we can use (9) for the skewness ratio at small ξ after replacing λ with λ + tα ′ v . For b = 2 this gives . (16) For small t we can write and thus approximate (16) by whereB Turning to the gluon distribution, we take with the function which has one parameter less than its counterpart (14).
Since most phenomenological information about gluons presently comes from small-x data, it would be difficult to constrain a third parameter. The analog of (16) reads for b = 2 and was already used in [28]. With the approximation in (17) we find wheref g (ξ) = α ′ g ln andB in analogy to the quark case. For our numerical study we take the parameters α ′ g = 0.164 GeV −2 and B g = 1.2 GeV −2 in order to match recent H1 data on J/Ψ photoproduction, whose t dependence is well fitted by with values b 0 = 4.63 GeV −2 and α ′ g = 0.164 GeV −2 for W 0 = 90 GeV [29]. To connect (27) with (24) we have used the approximate relation dσ/dt ∝ |H g (ξ, ξ, t)| 2 , which is obtained at tree level when one keeps only the imaginary part of the scattering amplitude. The skewness variable is given by 2ξ = (M J/Ψ /W γp ) 2 in terms of the γp c.m. energy. For simplicity, we have omitted the terms with 1/(n − λ) in (26) when fixing B g . For typical values of λ they are quite small.

For antiquarks we set
with x > 0. Little is known to date about the t dependence in the sea quark sector. Constraints can be provided by deeply virtual Compton scattering [30,31], which at small x is sensitive to both sea quark and gluon distributions. A comprehensive analysis of this data, as has recently been performed in [20], is beyond the scope of this work. We will instead explore the pattern of evolution for two extreme choices. In model 1 we set the t slope fq equal to the one for valence quarks: where the choice fs = f d has no strong motivation, but does not strongly influence the results we will obtain. In model 2 we set instead for all quark flavors. The initial conditions for evolution of the singlet and non-singlet combinations are then obtained from For the evolution study in the next section, we make the ansatz (3) with the CTEQ6M parton distributions at µ 0 = 2 GeV, so that we can use the fit of [27] for the t dependence of H q v (x, 0, t) as specified at the beginning of this section. In (31) we have neglected the tiny charm quark distribution at µ 0 . To explore the region of lower scales, we will also consider backward evolution.

Evolution of the t dependence
In accordance with the analytical considerations in the previous section, we find that at the initial scale the t dependence of H g (ξ, ξ, t) is well described by an exponential form at small t and ξ. Evolving to higher scales we still find an approximately exponential behavior in both model 1 and 2, as shown in Fig. 9 for ξ = 3.2 × 10 −4 . A slight departure from an exact exponential in the full region 0 ≤ −t ≤ 1 GeV 2 is however visible at µ 2 = 50 GeV 2 . Evolving to lower scales, we still find an approximate exponential t dependence at µ 2 = 3 GeV 2 , but for yet lower scales the situation changes. At µ 2 = 2 GeV 2 the distribution H g (ξ, ξ, t) turns negative for −t around 0.5 GeV 2 in model 1 and around 0.3 GeV 2 in model 2, whereas at µ 2 = 1.69 GeV 2 we have H g (ξ, ξ, t) < 0 already for t = 0. This is due to the behavior of the CTEQ6M gluon density at low scales. Since the gluon distribution in this region varies considerably between different global parton fits, we shall not elaborate on this issue further here.
The singlet distribution H S (ξ, ξ, t) is again well approximated by an exponential in t at the starting scale, and it stays exponential to high accuracy in model 2 up to µ 2 = 50 GeV 2 and even down to µ 2 = 2 GeV 2 . As shown in Fig. 9, this is however not the case in model 1. Here we find a clear departure from an exponential behavior even when evolving from the starting scale to µ 2 = 6 GeV 2 , whereas under backward evolution H S (ξ, ξ, t) rapidly turns negative for some value of t. We notice that in model 1 the x dependence of H S (x, 0, t) at the starting scale rapidly changes with t due to the large value of α ′ v . This induces a corresponding change in the x dependence of H S (x, ξ, t), which enters in the evolution equations.
To quantify the change of the t dependence under evolution, we fit the GPDs at given ξ and µ to for −t between 0 and 0.5 GeV 2 , where i = g, S. Given the behavior of the distributions under backward evolution, we restrict these fits to µ 2 ≥ 4 GeV 2 . Whereas for H S in model 2 and for H g in both models the form (32) gives an excellent description in the kinematical region of the fit, the corresponding fit for H S in model 1 can only be approximate, as is seen in Fig. 9. This must be kept in mind when interpreting the subsequent results, but despite this caveat the corresponding t slopef S gives a fair account of how H S (ξ, ξ, t) changes with µ. The results of the fit are shown in Fig. 10 for the starting scale and for µ 2 = 50 GeV 2 . We see that over a wide region of small ξ the dependence off i (ξ; µ) on ξ remains logarithmic after evolution to higher scales. For given µ we can hence perform a fit The results of such a fit in the range 3.2 × 10 −5 < ξ < 3.2 × 10 −4 are shown in Fig. 11, where we plot f i (ξ; µ) at the midpoint ξ = 10 −4 of the fit range, as well as the effective shrinkage parameter α ′ i (µ). In model 2,f i (ξ; µ) and α ′ i (µ) are equal for the gluon and the quark singlet to a good precision at the starting scale by construction. They change rather mildly under evolution to higher scales, but a visible difference between gluon and singlet appears, especially for α ′ i .   Figure 10: The t slopef (ξ; µ) fitted according to (32) for the gluon (left) and the quark singlet GPD (right). At µ 2 = 4 GeV 2 the curves forf g (ξ; µ) coincide in models 1 and 2. In model 1, we see that the slope and the shrinkage parameter for the singlet evolve quite strongly and tend to approach the corresponding values in the gluon distribution, which increasingly dominates evolution with increasing µ. The respective values for the gluon and the quark singlet are however clearly different even at µ 2 = 50 GeV 2 .
We have also investigated the evolution behavior of the non-singlet quantity H u−d (ξ, ξ, t). Only the flavor difference enters for the sea quark distributions (28) in this case, and we restrict our investigation to model 1. The t dependence changes in a similar way as for the quark singlet in model 1, as becomes evident from comparison of Fig. 9 and Fig. 12. In particular, evolution to higher scales modifies the exponential behavior of the initial condition. A fit of the t slope for −t between 0 and 0.5 GeV 2 must hence be taken with the same caveat as above. The t slopef u−d (ξ; µ) fitted as in (32) shows an approximately logarithmic ξ behavior in the full range 2 GeV 2 ≤ µ 2 ≤ 50 GeV 2 , so that we can again perform a fit to the form (33) in the region 3.2 × 10 −5 < ξ < 3.2 × 10 −4 . Between µ 2 = 2 GeV 2 and 50 GeV 2 , the resulting shrinkage parameter α ′ u−d increases only by about 2%, and the slopef u−d (ξ; µ) at ξ = 10 −4 decreases by about 10%. Compared with the quark singlet sector in model 1, evolution effects in the non-singlet sector are hence considerably weaker.

Conclusions
We have studied several aspects of the evolution behavior of GPDs at small ξ. To do this, we assumed a particular form of the GPDs at a moderately low scale and numerically evolved this model ansatz to higher µ. At t = 0 we have taken initial conditions for which H(ξ, ξ, 0) and the corresponding parton den-sity H(ξ, 0, 0) approximately obey power-laws with the same power. Under evolution to higher scales this power changes but remains the same for a GPD and the associated PDF. As a consequence, the skewness ratio R(ξ, µ) is only weakly ξ dependent, to the extent that the effective power changes with ξ. The values of R(ξ, µ) for different initial conditions approach each other with increasing µ, and at high scales they are well approximated by the Shuvaev formula (10). This convergence is however not very fast: with rather different values of R(ξ, µ) at µ 2 = m 2 c it only becomes visible at µ 2 of a few 10 GeV 2 for the gluon and quark singlet distributions, and at yet larger values in the non-singlet sector. We have not attempted to study how the situation would change for initial conditions at much lower scale, considering that in this case the leading-order approximation of the evolution equations would no longer be suitable for drawing quantitative conclusions. We confirm the finding of [10] that evolution to higher scales generates a singular derivative (∂/∂x)H(x, ξ, t) at x = ξ for quarks, but not for gluons.
To study the change of t dependence under evolution, we have chosen initial conditions at µ 0 = 2 GeV such that H(ξ, ξ, t) ∼ exp[tf (ξ)] andf (ξ) = α ′ ln(1/ξ) +B at small ξ and t. For distributions with a small shrinkage parameter α ′ , we find that to a good approximation the t dependence remains exponential under evolution to higher (and to some extent also to lower) scales. In contrast, a deviation from an exponential t behavior becomes visible after evolution rather quickly for distributions with large α ′ at the starting scale, so that a fit to an exponential form is only approximate in these cases. The fitted slopes f (ξ) of the evolved GPDs retain a logarithmic ξ dependence, so that one can also determine a shrinkage parameter α ′ at different scales µ. We find that the values of α ′ for the gluon and quark singlet distributions remain close (but not equal) to each other under evolution in a model where they coincide at µ 0 . In an alternative model, where α ′ for the quark singlet is much larger than for gluons at µ 0 , evolution brings their values closer to each other, but clear differences remain even at µ 2 = 50 GeV 2 . An analogous behavior is found forf (ξ) at given ξ in both models. We therefore conclude that one may not take it for granted that the t dependence of gluon and sea quark distributions is the same at moderate scales. In the flavor nonsinglet sector, we find thatf (ξ) and α ′ remain quite stable under evolution of the scale.