Kinematical higher-twist corrections in $\gamma^* \gamma \to M \bar M $

We estimate kinematical higher-twist (up to twist 4) corrections to the $\gamma^*(q_1) \gamma(q_2) \to M(p_1) \bar{M}(p_2)$ amplitudes at large $Q^2=-q_1^2$ and small $s=(q_1+q_2)^2$, where $M$ is a scalar or pseudoscalar meson. This process is known to factorize at leading twist into a perturbatively calculable coefficient function and generalized distribution amplitudes (GDAs). The kinematical higher-twist contributions of order $s/Q^2$ and $m^2/Q^2$ turn out to be important in the cross section, considering the kinematics accessible at Belle and Belle II. We present numerical estimates for the cross section for $\gamma^* \gamma \to \pi^0 \pi^0$ with the $\pi \pi$ GDA extracted from Belle measurements and with the asymptotic $\pi \pi$ GDA as inputs to study the magnitude of the kinematical corrections. To see how the target mass corrections of order $m^2/Q^2$ affect the cross section, we also perform the calculation for $\gamma^* \gamma \to \eta \eta$ by using a model $\eta \eta$ GDA.In the range $s>1$ GeV$^2$, the kinematical higher-twist corrections account for $\sim 15 \%$ of the total cross section, an effect which is not negligible. Since $\pi \pi$ GDAs are the best way to access the pion energy-momentum tensor (EMT), our study demonstrates that an accurate evaluation of EMT form factors requires the inclusion of kinematical higher-twist contributions.

As the studies of GPDs allow us to perform nucleon tomography through a Fourier transform in the transverse coordinate space [14][15][16], GDAs open the way to an impact-parameter picture [17] of the exclusive hadronization process q q → M (p 1 ) M (p 2 ).The GPDs and GDAs are also used to investigate the matrix elements of the energymomentum tensor (EMT) [3,[18][19][20] for hadrons in the spacelike and timelike regions, respectively.One can extract mass, pressure and shear force distributions of hadrons with the spacelike EMT form factors [21][22][23][24][25][26][27][28].Since there is no experimental facility where pion GPDs can be directly measured (see however Refs.[29][30][31]), the studies of ππ GDAs are a necessary tool to access the pion EMT.The spacelike pion EMT form factors can be obtained from the timelike ones by using dispersion relations, and in this process the ππ GDAs and pion EMT form factors at s > 1 GeV 2 should be included so as to make the integrals convergent.Therefore, this goal necessitates to extract GDAs in a sufficiently large s-range, thus demanding a control as precise as possible of kinematical higher-twist corrections to the amplitudes, which are proportional to s/Q 2 and m 2 /Q 2 , with m the meson mass.While the contribution of one higher-twist process has been previously discussed in Ref. [32], it only matters in a very limited kinematical region, namely the near forward or near backward regions.The phenomenological necessity of a sizeable (genuine) twist-4 contribution to the γ * γ → ρρ amplitude has also been pointed out [33].A complete understanding of higher-twist corrections to the process γ * (q 1 )γ(q 2 ) → M (p 1 ) M (p 2 ) , is however a difficult task which is far from being achieved.Meanwhile, a separation of kinematical and dynamical contributions in the product of two electromagnetic currents T {j em µ (z 1 x)j em ν (z 2 x)} was proven in Refs.[34][35][36] and applied to the deeply-virtual Compton scattering (DVCS) reaction [37].The kinematical corrections come from two types of operators, namely the subtraction of traces in the leading-twist operators and the higher-twist operators which can be reduced to the total derivatives of the leadingtwist ones.The subtraction of traces was applied in Ref. [38] to the reaction of Deep Inelastic Scattering (DIS), leading to target mass corrections.The kinematical corrections in DVCS can be considered as a generalization of these target mass corrections.However, higher-twist operators which can be reduced to the total derivatives of the leading-twist ones will also contribute to the DVCS reaction, since nonforward matrix elements are used.As pointed out in Refs.[34][35][36], the distinction between two types of kinematical corrections is not Lorentz invariant and has no physical meaning.Both contributions should therefore better always be added together.Since the same operator governs the physics of the reaction (1), one may use the same techniques to improve our understanding of the sdependence of its amplitude.We thus study here the kinematical higher-twist corrections to the amplitude of the reaction (1), in the kinematical domain suitable for a collinear QCD factorization framework where the leading-twist amplitude can be written as the convolution of a perturbatively calculable coefficient function and GDAs [20].
In Sect.II, we describe the kinematics of the γ * γ → M M process and recall the basic properties of GDAs.In Sect.III, we recall the results of Refs.[34][35][36] and the definitions of the higher-twist kinematical operators.In Sect.IV, we derive the helicity amplitudes for the reaction (1), including the kinematical higher-twist contributions.Sect.V shows our numerical estimates of the kinematical higher-twist contributions to the cross section for both ππ and ηη cases.We briefly present our conclusions in Sect.VI.Appendices A and B provide technical details for the calculation of helicity amplitudes.

II. KINEMATICS AND GENERALIZED DISTRIBUTION AMPLITUDES
FIG. 1: Kinematics of the process γ * (q1)γ(q2) → M (p1) M (p2) in the center of mass of the meson pair; the virtual photon is emitted by the electron with four-momentum q1 = k1 − k2.
To describe the process (1), we define the lightlike vectors n and ñ in a convenient way so that they can be expressed by the momenta of the spacelike virtual photon q 1 and the real photon where The polar angle of the meson (M ) momenta θ is illustrated in Fig. 1, and is defined as where m is the meson mass.For convenience, a new variable but the final amplitudes will be expressed in terms of cos θ.If the z-axis is chosen so that the t-z plane contains the lightlike vectors n and ñ, then only ∆ = p 2 − p 1 has a transverse momentum, ∆ = ζ 0 (ñ − τ n) + ∆ T .Using the on-shell condition, we obtain ∆ 2 The amplitude for γ * γ → M M is defined as where r = z 1 q 1 + z 2 q 2 , and the constraint z 1 − z 2 = 1 is imposed for real constants z 1 and z 2 .Owing to the electromagnetic gauge invariance, one can decompose this amplitude as [37] with g µν ⊥ and ǫ µν ⊥ given by The last term in Eq. ( 6) is of no interest since it does not contribute to any observable, and the rest of them can be expressed in terms of the GDAs if the factorization conditions Q 2 ≫ s, Λ 2 QCD are satisfied.The leading-twist amplitude was first presented in Ref. [20] with the help of a twist-2 GDA Φ q (z, ζ 0 , s) for an isoscalar meson pair, where P = (p 1 + p 2 )/2, Φ q is the GDA for the quark flavor q, q(z 1 n)/ nq(z 2 n) is the leading-twist vector operator (a light-like Wilson line joining the points z 1 n and z 2 n is implied), and z 1 − z 2 = 1 is not a necessary condition.This matrix element can alternatively be expressed in terms of double distributions (DDs) as [39] M (p 2 )M (p with f q and g q having support on the rhombus |α| + |β| ≤ 1 and assumed to vanish at the boundary, and Then, one can easily relate the GDA to double distributions where y = 2z − 1.Since the meson pair is produced with charge conjugation C = +1, one can obtain the relations from charge conjugation invariance.Assuming that the DDs vanish at the boundaries, Eq. ( 9) can be put in the form where the notation is introduced with symmetry φ q (β, α) = φ q (β, −α) = φ q (−β, −α), in order to simplify the calculation of the amplitudes thanks to the property where a, b and c are constants which are independent of α and β, and the exponents n and m are odd numbers.
Although the intermediate calculations involve the DD φ q (α, β), the final results will be presented in terms of the GDA using

III. OPERATOR PRODUCT EXPANSION AND HELICITY AMPLITUDES
A separation of kinematical and dynamical contributions in the time-ordered product of two electromagnetic currents i T {j em µ (z 1 x)j em ν (z 2 x)} was recently proved in Refs.[34][35][36].The kinematical contributions only involve the leadingtwist distributions, whereas unrelated genuine higher-twist distributions are necessary for the dynamical contributions.One can thus improve the description of reactions where two photons are involved by including the kinematical corrections, without any knowledge of the higher-twist distributions.A complete calculation of kinematical corrections was performed up to the twist-4 accuracy for DVCS with a (pseudo)scalar target in Ref. [37].In this work we shall apply similar techniques to calculate the kinematical higher-twist contributions in the s-t crossed channel of DVCS, namely the reaction γ * γ → M M .The kinematical contributions to the operator i T {j em µ (z 1 x)j em ν (z 2 x)} were given to twist-4 accuracy by [34,35,37], where the convention ǫ 0123 = 1 is adopted for the antisymmetric tensor and S µανβ is defined as In Eq. ( 17), V µ and W µ contain contributions of twist 2, twist 3 and twist 4, whereas X and Y are purely twist 4, see Appendix A for the detailed expressions.In practice, the spinor formalism [40,41] is used to calculate the amplitudes, since the expression of T µν becomes more compact and it is easier to figure out the twist of each term in the corresponding matrix elements.
In order to calculate the helicity amplitudes of Eq. ( 6), the photon polarization vectors are required.Choosing the momentum of the virtual photon along the z-axis, as shown in Fig. 1, its polarization vectors read [20] where the lower indices ± and 0 indicate the helicities of the photon.The polarization vectors ǫ of the real photon only have the transverse components, and they are related to the ones of the virtual photon as ǫ± = −ǫ ∓ .In the reaction γ * γ → M M , the helicity amplitudes are defined as and there are only three independent helicity amplitudes owing to parity invariance, as one can check from Eq. ( 6).
Here we choose the independent helicity amplitudes as A ++ , A 0+ and A −+ , then one obtains At leading twist, the operator product expansion of i T {j em µ (z 1 x)j em ν (z 2 x)} leads to the nonlocal operator with a lightlike separation.Since we are interested in the reaction γ * γ → M M with a charge conjugation even final state, one can safely neglect the contribution of the strange quark in the case of a π meson pair, where χ = 5e2 /18 is obtained thanks to the isospin symmetry1 .One needs however to add e 2 s s(z 1 n)/ ns(z 2 n) to the operator O ++ (z 1 n, z 2 n) in the case of a K meson pair.The kinematical higher-twist contributions in the operator product expansion of i T {j em µ (z 1 x)j em ν (z 2 x)} are related to the operator O t=2 ++ (z 1 , z 2 ), where the separation x is now not necessarily lightlike.We thus need to use the leading-twist projector Π(x, n) defined in Refs.[34][35][36], Since the dependence on n is always carried by a function of the type e −il•n in Eq. ( 13), the action of the leading-twist projector is simply given by Up to 1/Q 2 -accuracy, one obtains where φ = φ u + φ d and the second term provides a twist-4 contribution.In addition to the leading-twist operator O t=2 ++ (z 1 x, z 2 x), there are also the higher-twist operators which contribute to kinematical higher-twist corrections.Using Eq. ( 26), the matrix elements of O 1 and O 2 can be expressed up to 1/Q 2 -accuracy as Since the operators O 1 and O 2 contain total derivatives, their matrix elements vanish in the forward limit and need not be considered in DIS.They provide however corrections of order m 2 /Q 2 and s/Q 2 in the reaction γ * γ → M M .

IV. HELICITY AMPLITUDES IN TERMS OF GDAS
In the following we will calculate the helicity amplitudes of γ * γ → M M , adopting similar techniques to the ones used for DVCS in Ref. [37].There are three independent helicity amplitudes, which can be expressed in terms of DDs, where ∆ 2 T = g µν ⊥ ∆ µ ∆ ν and iǫ is omitted in the functions of ln and Li 2 since it will not contribute to the amplitudes.Details of the calculations can be found in Appendix B. A 0+ and A −+ are proportional to ∆ • ǫ − and (∆ • ǫ − ) 2 as indicated by Eq. ( 21), respectively, and the amplitudes do not depend on z 1 and z 2 which indicates that the translation invariance is recovered in the physical amplitudes.The function F (α, β) is defined as where F = 0 and F = 1 correspond to the quark momentum fractions z = 1 and z = 0 of the GDAs, respectively.We notice that there are three types of integrals expressed by DDs in the obtained amplitudes, namely where Y (F ) is some function of F .Inserting the identity dy δ(βζ 0 − y − α) = 1 into the integrals above, one can reexpress the integrals in terms of GDAs by using Eq. ( 16), where y = 2z − 1 and ζ 0 = −β 0 cos θ as defined in Eq. ( 3).Therefore, we can write the helicity amplitudes as where η = cos θ and Φ = Φ u + Φ d .The GDA for s quarks is also required in some reactions such as γ * γ → K 0 K0 with a charge conjugation-even K meson pair, and we just need to replace χΦ with e 2 u Φ u + e 2 d Φ d + e 2 s Φ s in the above amplitudes.One can clearly see the O(s/Q 2 ) corrections in the amplitudes, and the target mass correction of order In general, charge conjugation-even GDAs can be expanded as [20] Φ(z, cos θ, s where C (3/2) n (x) are Gegenbauer polynomials and P l (x) are Legendre polynomials.Due to this general expression for GDAs, the singularities of 1 z , 1 1−z and ln (1 − z) in the helicity amplitudes will be compensated by the GDA when z → 0 and z → 1.As a consequence, the amplitudes have no end-point singularities.In the asymptotic limit (Q 2 → ∞), only the terms with n = 1 survive, where the first and second terms correspond to the S-wave and D-wave production of a meson pair, respectively.The nonvanishing helicity-flip amplitudes A −+ (A (2) ) and A 0+ (A (1) ) indicate the existence of a D-wave GDA.

V. NUMERICAL ESTIMATES OF THE HIGHER-TWIST KINEMATICAL CONTRIBUTIONS
A. ππ GDA extracted from Belle measurements The process γ * γ → M M can be measured in e + e − collisions, which are accessible at KEKB and SuperKEKB.In Ref. [20], the differential cross section for eγ → eM M is expressed as where ϕ is the azimuthal angle of the meson pair as illustrated in Fig. 1 and s eγ is the center-of-mass squared energy of eγ.ǫ is defined as usual by 2: Differential cross section for eγ → eπ 0 π 0 calculated with ππ GDA extracted from Belle measurements through a leading-twist analysis [9].Dashed curves show the twist-2 results, while solid curves include the kinematical higher-twist contributions.The selected values are seγ = 30 GeV 2 , Q 2 = 9 (16, 25) GeV 2 and cos θ = 0.2 (0.4, 0.6, 0.8) as indicated on the different panels.
In the reaction eγ → eππ, there are two types of contributions to the cross section.The final π + π − with negative charge conjugation couples to a virtual photon, and its contribution is expressed in terms of the pion electromagnetic form factor.When the charge conjugation of ππ is positive, the pion pair can be π 0 π 0 or π + π − .Using factorization, this type of contribution is determined by GDAs, which we are interested in.The π + π − GDA is equal to the one of π 0 π 0 due to the isospin symmetry.However, since π 0 π 0 are identical bosons, cos θ will be restricted to 0 ≤ cos θ ≤ 1 in Eq. (36).After integration over θ and ϕ, the cross section for C-even π + π − production is then twice as large as the one for π 0 π 0 .
In 2016, the Belle Collaboration released the measurements of differential cross section for γ * + γ → π 0 + π 0 [8].Since the final state is π 0 π 0 , there is no contribution from the pion electromagnetic form factor.The twist-2 ππ GDA was extracted by using the leading-twist amplitude [9].We use this pion GDA to estimate the cross section for eγ → eππ where the integral over ϕ is performed in Eq. (36), In order to show the size of the higher-twist kinematical contributions, Eqs. ( 33) and ( 38) are used to calculate the cross section, and the results are depicted as the solid lines in Fig. 2. The dashed lines represent the leading-twist cross sections.Considering the kinematics of Belle measurements, we choose the values Q 2 = 9, 16, 25 GeV 2 , s ∈ (0.25, 4) GeV 2 , and we set s eγ = 30 GeV 2 which is the typical value at Belle.In Fig. 2, black lines denote cos θ = 0.2 and orange lines correspond to cos θ = 0.4, while cos θ = 0.6 and cos θ = 0.8 are depicted as red and blue, respectively.As Q 2 increases, kinematical contributions become less important, which is consistent with the fact that the kinematical contributions are suppressed by 1/Q or 1/Q 2 .The kinematical contributions cannot be neglected in the region where √ s ≥ 1 GeV.The helicity-flip amplitudes A −+ and A 0+ receive only contributions from the D-wave GDA, and a large difference between two types of cross sections is displayed around the D-wave resonance region of f 2 (1270) in Fig. 2. Hence, the study of the amplitudes A −+ and A 0+ will be important for the investigation of this resonance region.The kinematical higher-twist corrections contribute ∼ 15% to the cross section on average, if one restricts the process eγ → eππ to the kinematics of Belle measurements.
In Fig. 3, we also present the ratio dσ(2 + 3 + 4)/dσ(2) where dσ(i) (i = 2, 3, 4) is the twist-i contribution to the cross section, and the colors of the lines indicate different values of cos θ as in Fig. 2. In this figure, the contributions of the kinematical higher-twist corrections are quite clear, so we can infer that the kinematical corrections cannot be neglected when √ s > 1 GeV.Around √ s ∼ 1.5 GeV, the kinematical corrections are dominant in the cross section with cos θ = 0.8; this appears because the twist-2 cross section is quite tiny when calculated with the GDA extracted from Belle measurements; this GDA may however not be accurate in this region since the uncertainties of Belle measurements are quite large there; this ratio may thus not reflect the real physics around √ s ∼ 1.5 GeV.As we have seen, the kinematical corrections are not negligible in the region √ s > 1 GeV, which turns out to be important for the studies of the pion EMT form factors. Indeed, since pion GPDs cannot easily be measured in experiments, GDAs offer a way to investigate the timelike EMT form factors of pions.The spacelike EMT form factors can then be obtained from the timelike ones by using dispersion relations, in which case the timelike EMT form factors of √ s > 1 GeV are needed to be included numerically.As a consequence, it is important to use the most accurate description of the cross section with the inclusion of kinematical contributions.
As pointed out above, the uncertainties of Belle measurements [8] are quite large, and the statistical errors are dominant.However, this situation will be improved substantially soon, since the Belle II collaboration just started taking data at the SuperKEKB with a much higher luminosity.Precise measurements of γ * +γ → M + M are expected in the near future, and an accurate description of the amplitudes for the study of GDAs requires the inclusion of kinematical contributions up to twist 4. The asymptotic pion GDA used in our calculation is taken from Eq. (68) of Ref. [2], where δ 0 and δ 2 are ππ elastic scattering phase shifts in the isospin 0 channel [42][43][44].R π = 0.5 represents the momentum fraction carried by quarks in the pion meson.In this asymptotic GDA, we do not include the contribution of the f 2 resonance.However, we believe it is reasonable to use this GDA here, since our purpose is not to predict the cross section for eγ → eπ 0 π 0 precisely, but to estimate the magnitude of kinematical higher-twist contributions and determine whether one can neglect them or not in the cross section.In Fig. 4, we show the cross section for eγ → eπ 0 π 0 with fixed Q 2 and cos θ, the dashed lines are the twist-2 cross sections, while the solid ones indicate the cross sections with kinematical contributions included.The colors of the lines denote different values of cos θ as indicated on the different panels of the figure.Similarly to the case of the extracted ππ GDA, the kinematical corrections are important to describe the cross section in the region of √ s > 1 GeV.As Q 2 increases, the kinematical contributions become less important.Compared with Fig. 2, the magnitude of the cross sections are similar at different Q 2 , even though the asymptotic GDA is very different from the extracted

VI. SUMMARY
In this paper we presented a complete calculation of kinematical higher-twist corrections for the helicity amplitudes of the reaction γ * + γ → M + M up to twist 4, where only the leading-twist GDA is involved in the description of the cross section.In case of ππ production, we use two types of GDAs to estimate the kinematical higer-twist contributions in the cross section, namely the ππ GDA extracted from the Belle measurements and the asymptotic ππ GDA.Even though those two GDAs are very different, both of them lead to kinematical corrections which cannot be neglected for γ * + γ → ππ considering the kinematics accessible at Belle and Belle II.Moreover, the relative magnitude of the higher-twist corrections is also comparable for the two types of ππ GDAs at different values of Q 2 , as seen from Figs.
Due to the small pion mass, only kinematical corrections of the type O(s/Q 2 ) contribute to the cross section for γ * + γ → ππ.The production of a pair of slightly heavier mesons is needed to check how the target mass corrections m 2 /Q 2 affect the cross section.Since the ηη GDA is an unknown quantity at present, we calculate this effect with a model GDA identical to the asymptotic ππ GDA except that the mass of η is used.This calculation indicates that kinematical corrections are also not negligible in this case.Furthermore, the negative kinematical corrections from the amplitude A ++ are dominant over the higher-twist contributions from A +− and A 0− , and the kinematical corrections are always negative in the cross section, which is different from the ππ case where kinematical corrections can go both ways.In conclusion, let us stress that while the uncertainties of present Belle measurements [8] are too large for our study to invalidate the conclusions of [9], the situation should change substantially in a near future since Belle II collaboration just started taking data at the SuperKEKB with a much higher luminosity.An accurate description of the amplitudes in terms of GDAs will therefore require the inclusion of kinematical contributions up to twist 4. Note also that precise measurements of γ * + γ → M + M for various mesons will be of utmost importance to address the questions of the pion EMT form factors and of the impact-parameter representation of GDAs [17].
In the future, this work can be extended to the production of other meson pairs, for example the γ * γ → πη channel which should help unraveling the quark and gluon structure of hybrid meson (J P C = 1 −+ ) [46].The scattering amplitude is likely to be sensitive to sizeable kinematical higher twist contributions.This channel has recently been advocated [47] to be related to shear viscosity of quarks in hadronic matter.The production of a pair of vector mesons should also be discussed, opening the way to a meaningful extraction of the EMT form factors for ρ or ω mesons.
Similar relations apply to the timelike process amplitude γ * → M M γ which opens another access to GDAs [48] through the interference with the initial-state radiation amplitude in the process e + e − → M M γ, as experimentally proven by the BABAR collaboration [49].The extension of our work to the process γ * + γ → N + N will also be needed if Belle II detector is able to detect this channel.
We consider the matrix element of the twist-2 part of T 0+ = T µν ǫ µ 0 ǫν + and substitute it into Eq.( 5), where (r + ul z1z2 ) 2 = −2 n • ñ (uF + (1 − u)z 1 ) is obtained by neglecting the terms of order O(s, m 2 ) since they will not contribute at the 1/Q 2 accuracy.Let us mention that iǫ can be omitted in ln(F − iǫ) since F is always positive (0 ≤ F ≤ 1) and there is no branch cut.Indeed, ln(F ) is divergent when F = 0 (z = 1), but the GDA vanishes at z = 1.Similarly, the twist-3 part is obtained, We note that although the twist-2 and twist-3 amplitudes depend on z 1 and z 2 , this dependence disappears in their sum, which indicates that translation invariance is recovered in the physical amplitudes.Once again iǫ is omitted in ln(F − iǫ) as well.The calculation of A −+ is quite similar to the one of A 0+ , For the sum, one gets which does not depend on z 1 and z 2 , and which is proportional to (∆ • ǫ − ) 2 as shown in Eq. ( 21).The calculation of A ++ is more lengthy than the ones of helicity-flip amplitudes, since it contains the contributions of twist 2, twist 3 and twist 4.
Taking the trace of Eq. ( 6), one obtains [37] A We simplify the operator expansion of T ++ in Eq. ( 17) by using Eq.(B9), and we take the matrix elements of T ++ to obtain the contribution of each twist to A ++ [37], where the expressions for B t=3 µ (z 1 , z 2 ) and R(z 2 , z 1 ) can be found in Appendix A. The calculation of A t=2 ++ is rather straightforward with the help of Eq. ( 26), where Eq. ( 15) is used to eliminate the irrelevant terms, and the first term in Eq. ( B11) is actually the twist-2 amplitude given by Ref. [2].The twist-3 contribution is expressed as and is independent of z 1 and z 2 .There are several terms contributing to the amplitude of A t=4 ++ in Eq. (B10), which we denote as 4) .(B13) A t=4 (1) is calculated by using Eq. ( 28), When calculating the other terms in A t=4 ++ , one finds divergences in A t=4 (2) and A t=4 (3) associated with the real (onshell) photon.Here, the photon is set as an offshell one so as to regularize the divergences, q 2 → q 2 + (ξ − z 1 )(q 1 + q 2 ), r = z 1 q 1 + z 2 q 2 → r = −q 2 + ξ(q 1 + q 2 ), (B15) and we take ξ = z 1 at the end of the calculations to obtain the final results.As shown in Eq. (B10), the operator R(z 1 , z 2 ) is involved in the remaining terms in A ++ t=4 .First, we simplify the matrix element of R(z 1 , z where z 12 is kept so that z 12 = 1 is not a necessary condition.One can make the replacements z 1 → uz 1 and z 2 → uz 2 to obtain R(z 1 u, z 2 u).We substitute the matrix elements R(z 1 u, z 2 u) and R(z 2 u, z 1 u) into Eq.(B10) and we take ξ = z 1 , Similarly, the last contribution to A t=4 ++ is given by Summing over all contributions leads finally to = χ dβ dα φ(β, α) 2 ln(F ) + s n • ñ (F − α) + FIG.2: Differential cross section for eγ → eπ 0 π 0 calculated with ππ GDA extracted from Belle measurements through a leading-twist analysis[9].Dashed curves show the twist-2 results, while solid curves include the kinematical higher-twist contributions.The selected values are seγ = 30 GeV 2 , Q 2 = 9 (16, 25) GeV 2 and cos θ = 0.2 (0.4, 0.6, 0.8) as indicated on the different panels.

FIG. 6 :
FIG.6: Differential cross section for eγ → eηη with the model ηη GDA described in the text, same conventions as in Fig.2.