Twist-3 contributions to $\gamma\gamma\rightarrow\pi^+\pi^-,K^+K^-$ processes in perturbative QCD approach

As one of the simplest hadronic processes, $\gamma\gamma\rightarrow M^{+}M^{-}$ ($M=\pi,K$) could be a good testing ground for our understanding of the perturbative and nonperturbative structure of QCD, and will be studied with high precision at BELLE-\RNum{2} in the near future. In this paper, we revisit these processes with twist-3 corrections in the perturbative QCD approach based on the $k_{T}$ factorization theorem, in which transverse degrees of freedom as well as resummation effects are taken into account. The influence of the distribution amplitudes on the cross sections are discussed in detail. Our work shows that not only the transverse momentum effects but also the twist-3 corrections play a significant role in the processes $\gamma\gamma\rightarrow M^{+}M^{-}$ in the intermediate energy region. Especially in the few GeV region, the twist-3 contributions become dominant in the cross sections. And it is noteworthy that both the twist-3 result of the $\pi^{+}\pi^{-}$ cross section and that of the $K^{+}K^{-}$ cross section agree well with the BELLE and ALEPH measurements. For the pion and kaon angular distributions, there still exist discrepancies between our results and the experimental measurements. Possible reasons for these discrepancies are discussed briefly.


I. INTRODUCTION
A meaningful and historic subject of the perturbative QCD, which has received considerable attention in the past few decades, is the study of the exclusive processes at large momentum transfer. The pioneer works in this area are performed by Efremov and Radyushkin [1] and, independently, Brodsky and Lepage [2]. They pointed out that the exclusive processes at large momentum transfer can be factorized into perturbatively calculable kernels and hadronic wave functions. Although there is general agreement that perturbative QCD is able to make reliable predictions for the exclusive processes in the large energy region [3], the applicability of perturbative QCD to these processes in the intermediate energy region has been developed in controversy. In the field of exclusive processes, two-photon collisions γγ → π + π − , K + K − [4][5][6][7][8] are the specially ones with the initial states simple and controllable and the strong interactions only in the final states. These characteristics make the processes γγ → π + π − , K + K − be a good testing ground for our understanding of the perturbative and nonperturbative structure of QCD.
In the perturbative QCD approach based on the collinear factorization, the first investigation of the processes γγ → π + π − , K + K − was carried out by Brodsky and Lepage [9]. At leading-twist leading-order level, they gave detailed calculations of the helicity amplitudes of these processes and presented the known relation between the differential cross section and the timelike electromagnetic form factor: Several years later, Nižić [10] performed the leading-twist next-to-leading-order perturbative QCD calculations for the these processes. In Refs. [10,11], it was found that both the leading-order and next-to-leading-order perturbative QCD predictions are almost an order of magnitude smaller than the experimental data [12][13][14][15][16] in the intermediate energy region. Motivated by the calculation of the pion electromagnetic form factor [17] in which the nonleading twist contributions become the dominant one in a few GeV region, Gorsky [18] investigated the two-photon process γγ → π + π − with the higher-twist corrections. It was pointed out that the two-parton twist-3 contributions are logarithmically larger than the leading-twist contributions in the intermediate energy region due to the chiral enhancement effects. And this conclusion had also been confirmed in Ref. [19] with the BHL prescription [20]. But in a few GeV region, the predicted cross sections with higher-twist corrections [18,19] are several times larger than their experimental measurements [12][13][14][15][16].
As it is known that the higher-power effects from the intrinsic transverse momentum play a crucial role in the pion electromagnetic form factor at the scale of few GeV [21,22], one may expect that the same situation can also be found in the two-photon processes γγ → π + π − , K + K − . Based on the k T factorization theorem [21,23] in which the transverse momentum dependence was retained, the leading-twist perturbative QCD predictions for the cross sections of these two-photon processes were obtained in Refs. [24][25][26][27][28][29]. From the analysis of the differential cross sections obtained in both perturbative QCD and QCD sum rule, the conclusion that the transition from nonperturbative to perturbative QCD in γγ → π + π − , K + K − occurs at Q ≈ 2 GeV was drawn in Refs. [25,26,29]. However, the numerical results in Refs. [27,28] show that the twist-2 cross sections of the γγ → π + π − , K + K − in the k T factorization are still much smaller than the experimental data in the intermediate energy region.
To consider the higher-twist effects for the processes γγ → π + π − , K + K − in k T factorization would be a natural next step, and seems necessary because the higher-twist effects and the higher-power effects from the transverse momentum become important simultaneously in intermediate energy region. Compared with the previous calculations in collinear factorization [18,19], the higher-twist calculations in k T factorization have several major improvements: Firstly, both the transverse momentum dependence of hadronic wave functions and that of hard kernels are included. The transverse momentum dependence of the hard kernel regularizes the end-point singularities of the internal propagators. While the transverse momentum dependence of the wave function makes a sizable suppression to the perturbative QCD contributions, and cannot be ignored as pointed out in Ref. [22], especially in the few GeV region. Secondly, the corrections from the Sudakov resummation [21] substantially suppress the nonperturbative contributions from both the soft end-point regions and the large-b regions. Moreover, the corrections from the threshold resummation [30] suppress the nonperturbative contributions from the end-point regions further and eliminate the end-point singularity without a artificial cut-off [18] or the BHL prescription [19] in the twist-3 calculations. Thirdly, the scales of the coupling constant α s are chosen to be momentum-fraction dependent in order to avoid large logarithms from higherorder perturbative QCD corrections. All of the aforementioned improvements make the perturbative calculation become more self-consistent, even for momentum transfer as low as a few GeV.
In this paper, we present a detailed twist-3 calculations for the two-photon processes γγ → π + π − , K + K − in the perturbative QCD approach based on the k T factorization theorem. At the twist-2 level, our results are consistent with the predictions given in Refs. [27,28]. While, with the twist-3 corrections, it is found that both the π + π − and the K + K − cross sections agree well with the experimental data [13][14][15]. This paper is organized as follows: In Sect. II, the related conventions to these processes are introduced at first, next the twist-2 and twist-3 light-cone wave functions are described briefly, then we perform the calculation of the hard kernels with twist-3 corrections. The numerical analysis and discussion are given in Sect. III. The last section is our summary and conclusion. The expression of the Sudakov function can be found in Appendix.

A. Kinematics and conventions
We consider the simplest hadronic processes, γ 1 (p 1 , ε λ1 where the incoming photons are real with momenta p 1 , p 2 and polarization vectors ε λ1 1 , ε λ2 2 with the photon helicities λ 1 , λ 2 = ±1, M denotes the out-going pseudoscalar mesons with momenta p 3 , p 4 . The contributions to the γγ → (qq) + (qq) amplitude arise from 20 tree-level Feynman diagrams, which can be grouped together into four basic types (a, b, c and d) based on the virtualities of the internal gluon. We depict the four basic Feynman diagrams in Fig. 1. The remaining diagrams arise from permutations of the gluon and quark lines.
For convenience, the amplitude is calculated in the center-of-momentum frame with the momenta of the final mesons along z-axis. θ is the scattering angle and Q denotes the two-photon center-of-mass energy. In the light-cone coordinate, a momentum A µ has the form Under these conventions, the momenta of photons and mesons read and the corresponding polarization vectors of the initial photons can be written as with ω = Q 2 , s = sin θ 2 , c = cos θ 2 and p = (2ωsc, 0) for abbreviation. The momenta of the quark and antiquark k 1 , k 2 in the M + meson and k 3 , k 4 in the M − meson as labeled in Fig. 1 can be written as . Here x (y) andx (ȳ) are the longitudinal momentum fractions of quark and antiquark, respectively, and the transverse momenta k ⊥1 , k ⊥2 are relate to the momenta of the two out-going mesons.
The total number of tree-level diagrams in groups a, b, c, and d are 6, 6, 4, and 4, respectively. The grouping is based on the virtualities of the internal gluon.

B. Light-cone wave functions
The Fierz identity is generally employed to factorize the fermion flow: Here q andq are the quark and antiquark field. The Dirac structure γ µ γ 5 in Eq. (5) contributes at the twist-2 level, while iγ 5 and σ µν γ 5 contribute at the twist-3 level.
With the light-cone expansion to twist-3 accuracy, the two-parton light-cone wave functions of π − can be expressed via the matrix element of the bilocal operator [31][32][33][34] where f π refers to the decay constant of pion and the mass parameter µ π is defined as where m u , m d are the current quark masses and m π denotes the pion meson mass. Ψ π π is the twist-2 wave function and Ψ p π , Ψ σ π are two-parton twist-3 wave functions. Due to the SU (2) isotopic symmetry, the wave functions of π + are identical with the wave functions of π − . Generally, the twist-3 contributions are power-suppressed. However, compared with the twist-2 contributions, the suppression factor µ π /Q in the two-parton twist-3 contributions is usually not small enough because of the chiral enhancement effects, i.e., through the relations between light quark masses in the chiral perturbation theory, µ π (2 GeV) = 2.43 GeV [35]. Especially, in the intermediate energy region, such as the BELLE experimental range 2.4 GeV < Q < 4.1 GeV [15], one can obtain the suppression factor µ π /Q ∼ O(1). That is the reason why, in some processes [18,19,[36][37][38][39][40], the two-parton twist-3 contributions could become even more important than the twist-2 contributions in a few GeV region.
Transforming the matrix element of the bilocal operator in Eq. (6) from the coordinate space into the momentum space, one can obtain the light-cone projection operator [34,41,42] (8) with Ψ σ π = ∂Ψ σ π (x, k ⊥ )/∂x. From the first principles in QCD, the exact transverse momentum dependence of the above wave functions is still unknown. For the convenience of calculations, we assume that Ψ π π , Ψ p π and Ψ σ π have the same transverse momentum dependence. Following the ansatz made in Refs. [22,43], the wave functions of the pion can be separated into two parts One part is the distribution amplitudes φ π π , φ p π and φ σ π . Using the conformal symmetry, the distribution amplitudes can be expanded in the Gegenbauer polynomials [33] with the normalization conditions 1 0 dxφ π π (x, µ) = 1, Here a π 2,4 are the Gegenbauer moments, the twist-3 parameters η 3π , ρ 2 π are defined as and the parameters f 3π , ω 3π are defined in terms of matrix elements of the local three-parton twist-3 operators [33] 0|ūσ µα γ 5 g s G να d|π − (p) = 2if 3π p µ p ν , The other part is the k ⊥ -dependence function Σ π (x, k ⊥ ). With the constraint from the π − → µ −ν µ decay [20], one can obtain i.e., the function Σ π (x, k ⊥ ) satisfy the normalization condition An appropriate candidate for the function Σ π (x, k ⊥ ) is usually assumed to be a simple Gaussian form [22,44] where the oscillator parameter β π is determined by the mean square transverse momentum The detailed discussions about the choice of the oscillator parameter β π are shown in Refs. [22,44]. In this work, we will take the value k 2 ⊥ 1/2 π = 0.35 GeV, which is compatible with the π 0 → γγ constraint [20]. For the K − case, the matrix element of the bilocal operator can be expressed as [45] where f K refers to the decay constant of kaon, and the chiral enhancement parameter is defined as µ K = m 2 K /(m u + m s ), m u and m s are the current quark masses and m K denotes the kaon meson mass. We should note that the wave functions of K + are given by [45]. Similar to the pion case, the wave functions of the kaon meson can also be written as with Here we have adopted the approximation with the normalization conditions Here the SU (3) symmetry breaking effects have already been taken into account. The parameters a K 1,2 are the first two nonvanishing Gegenbauer moments, and the twist-3 parameters η 3K , ρ K ± read where the parameters f 3K , ω 3K , λ 3K are defined in terms of matrix elements of the local three-parton twist-3 operators [33] 0|ūσ µα γ 5 g s G να s|K − (p) = 2if 3K p µ p ν , Using the leading-order renormalization group equations, the scale dependence of the various parameters is given by for the pion and for the kaon with the evolution factor L = α s (µ 2 )/α s (µ 2 0 ) and β 0 = (33 − 2n f )/12. To the leading logarithmic accuracy, the QCD running coupling constant can be expressed as and the anomalous dimensions are given by with C F = 4/3 and C A = 3.
In perturbative QCD approach, the convolution of wave functions and hard scattering kernels needs to be performed in the transverse configuration b-space. With the definition one can obtain that the full soft wave functions of the pion and kaon in the transverse configuration b-space have the formŝ where the Fourier transformed functionsΣ M (x, b) read C. The helicity amplitudes with twist-3 contributions In the perturbative QCD approach, the helicity amplitudes can be expressed as the convolution between the soft hadron wave functions and the hard scattering kernels with respect to both the longitudinal momentum fractions and the transverse separation of quark and antiquark. To the π + π − production in two-photon collision, the helicity amplitudes have the form [26,29] T λ1λ2 ij represent the Fourier transformed hard scattering kernelŝ where T λ1λ2 ij with the transverse momentum k ⊥1 , k ⊥2 kept can be calculated perturbatively by Feynman rules. In the helicity amplitudes, the scripts i, j = π correspond to the twist-2×twist-2 term, and the scripts i, j = p, σ, σ correspond to the twoparton twist-3×twist-3 terms. It is noteworthy that the two-parton twist-2×twist-3 terms vanish because of the spin structure of the corresponding hard scattering kernels. Compared with the twist-2×twist-2 contributions, the contributions from the twoparton twist-3×twist-3 terms are generally regarded as power-suppressed. However in a few GeV region, the suppression factor is µ 2 π /Q 2 ∼ O(1) as a consequence of the chiral enhancement effects. So the two-parton twist-3×twist-3 contributions become important in this energy region. While, compared with the twist-2×twist-2 contributions, the contributions from the three-parton twist-3×twist-3 terms and that from the two-parton twist-2×twist-4 terms are suppressed by the factors f 3π /(f π Q) and m 2 π /Q 2 respectively. According to the values of the factors f 3π /(f π Q) and m 2 π /Q 2 in the order 10 −2 in the intermediate energy region, these two contributions are at least an order of magnitude smaller than the twist-2×twist-2 contributions in this energy region. A similar situation can also be found in the calculations of the pion and kaon electromagnetic form factors [38,40]. In this work, we will neglect the contributions form the three-parton twist-3×twist-3 terms and the two-parton twist-2×twist-4 terms.
In the following calculation, we will show our procedure to the Fourier transformation of the hard scattering kernels by means of a specific Feynman diagram in the Group a in Fig. 1 as an example. The hard scattering kernels T λ1λ2 a1ij with the subscript "a 1 " specifically refer to the results of one of the six diagrams in Group a. Its leading-twist contributions are illustrated as with the factor κ 1 = αα s (µ R )π 2 f 2 π . The corresponding quark and gluon propagators arẽ with The convolution formula in Eq. (32) involves two independent scales of this process, the factorization scale µ F and the renormalization scale µ R . In analogy to the cases of the pion electromagnetic form factor [21,37,46] and transition form factor [47][48][49], the transverse separation b 1 and b 2 of the pion wave functions provide the factorization scales µ F i = 1/b i (i = 1, 2), below which the QCD dynamics is regarded as being nonperturbative and can be absorbed into the the soft wave function. While the renormalization scale µ R = t is set to the largest mass scale associated with the momentum of the internal hard gluon propagators in order to minimize the higher-order QCD corrections.
To simplify the calculation, the hierarchy Q 2 xQ 2 ∼ yQ 2 xyQ 2 , k 2 1⊥ , k 2 2⊥ is postulated in the small-x, y region, as elaborated in Refs. [46,50]. Under this hierarchy, one can ignore the transverse momentum dependence of the internal quark propagators and retain it in the gluon propagator to regulate the endpoint singularity, then the hard kernels in Eq. (35) are simplified as where the terms contained the transverse momentum k ⊥ in the numerators, which are power-suppressed compared with other O(Q 2 ) terms [30], have also been dropped. By using the Fourier transformation of the propagators regularized with the i prescription one can obtain the Fourier transformed hard scattering kernels (b = |b 1 |) Here H 0 and K 0 denote Hankel and modified Bessel function. As a consequence of the above simplification, the hard kernels in Eq. (40) only depend on a single b parameter. The underlying physics picture is that, the virtual quark lines involved in the hard kernels are thought of as being far from mass shell, and shrunk to a point [21,26].
There are two types of resummation of the higher-order effects: the region with large transverse separation b for Sudakov resummation and the small longitudinal momentum fraction x region for threshold resummation. The Sudakov resummation of the double logarithms α s ln 2 [x/(Q 2 b 2 )] produced by overlapping collinear and soft divergencies for massless quarks as well as the renormalization group equation transformation from the factorization scale 1/b to the renormalization scale t, are incorporated into the Sudakov factor exp[−S]. In next-to-leading logarithm approximation, the Sudakov exponent S reads [21,26,51] S(x, y, b, where the function s(x, b, Q) is originally derived by Botts and Sterman [23] and later on slightly improved, for instance, by Li [21] and Kroll [51] et al. The explicit expression of s(x, b, Q) is given in Appendix. One can find that the Sudakov factor exhibits a strong fall off at large b. And this property makes the nonperturbative contributions from large b, no matter what x is, less important.
Since there exists the end-point enhancement in the two-parton twist-3 contributions (logarithmical enhancement) for the scattering process γγ → π + π − [18,19], the Sudakov factor exp[−S] is still not effective enough to suppress the nonperturbative contributions from small x region, which would spoil the perturbative calculation. It had been argued in Ref. [30] that as the end-point region is important, the corresponding large double logarithms α s ln 2 x that arise from the higher-order corrections also need to be organized to all orders, and into a jet function S t (x) as a consequence of threshold resummation. In next-to-leading logarithm accuracy, the jet function S t (x) can be parameterized into a universal form [30,52] In Ref. [48], Li and Mishima have proposed a parabolic parametrization for the parameter c, and frozen c at c = 1, when it exceeds unity: Since the factor S t (x) drops rapidly as x → 0, 1, the end-point singularities are eliminated and the nonperturbative contributions from the dangerous end-point regions are suppressed effectively. After including the above two types of resummation, both the twist-2 and twist-3 contributions are well-behaved, and that results in the perturbative calculation more self-consistent. Performing the integration over b 2 , the helicity amplitudes in Eq. (32) can be then rewritten as a much clearer and simpler form in a single-b form Making the substitution d 2 b → bdbdϕ and performing the integral over the polar angle ϕ, we can obtain the helicity amplitudes for γγ → π + π − process in the conjugate b space by summing up all the 20 Feynman diagrams with the twist-2 hard kernelsT T +− aππ = 16κ 1 81 4 + s 2 (xy + xȳ) c 2 xȳ + c 2 (xy + xȳ) s 2 xȳ F(s a , b), and the nonzero twist-3 hard kernelŝ Here the factor κ 2 is defined as Compared with the twist-2 hard kernels, the twist-3 ones are easily found to be suppressed by the factor κ 2 /κ 1 = µ 2 π /ω 2 . In the above expressions for theT λ1λ2 nij , J 0 represents the Bessel function with the variable pb = 2ωscb and the function F(s i , b) is defined as The renormalization scales t n (n = a, b, c, d), enteringT H as the arguments of strong coupling constant α s and depending on the kinematics of the specific group, are chosen as The factors s n (n = a, b, c, d) come from the gluon propagators in each group In addition, it is worth noticing the hard kernels with different helicity have the relationŝ For the γγ → K + K − process, the hard kernels can be obtained directly with the replacements of e d → e s , f π → f K , µ π → µ K .
And the corresponding helicity amplitudes can be expressed as (64)

III. NUMERICAL ANALYSIS
The differential cross sections of the processes γγ → M + M − (M = π, K) can be expressed as In order to make comparison with the experimental measurements, we integrate over the scattering angle in the region | cos θ| < 0.6 to get the cross sections σ 0 for both π + π − and K + K − processes. In this work, we focus on the intermediate energy region 1.0 GeV < Q < 7.0 GeV to check the applicability of the perturbative QCD in γγ → π + π − , K + K − processes. 0.25 ± 0.15 [35,53] π µ = 2 GeV K µ = 2 GeV units/Refs.
In one-loop accuracy, the QCD running coupling is given by α s (µ 2 ) = π/[β 0 ln(µ 2 /Λ 2 QCD )] with the QCD scale Λ QCD = 0.2 GeV and the interaction scale µ = t n (n = a, b, c, d) corresponding to each group in Fig. 1. The most important parameters used in our numerical analysis are listed in Table. I. In addition, the decay constants are taken as f π = 130.4 MeV and f K = 159.8 MeV. For the meson masses and the quark mass entering the pion and kaon distribution amplitudes, we take m π = 139.6 MeV, m K = 493.7 MeV and m s = 95 MeV in MS scheme. All of the values are quoted from PDG [55]. The Gegenbauer moments a π 2,4 and a K 1,2 in the pion and kaon distribution amplitudes come from QCD sum rule [35,53]. For the chiral enhancement parameters µ π and µ K , we employ the well known chiral perturbative theory relations [56] and obtain The remaining parameters of the twist-3 distribution amplitudes in our calculations are taken from Ref. [45] and shown in Table. I for brevity.  The theoretical uncertainty in this work mainly comes from the nonperturbative inputs of meson distribution amplitudes: one is the Gegenbauer moments a M n and the other is the chiral enhancement parameters µ M (n = 2, 4 for pion and n = 1, 2 for kaon, M = π, K). In Table II, we present five different models of the distribution amplitudes. And the corresponding Gegenbauer moments are in the region which can cover the ones given in the recent lattice determinations [57].
With the five models of the pion and kaon distribution amplitudes, the cross sections σ 0 (M + M − ) at twist-2 and twist-3 levels are illustrated in Fig. 2 and Fig. 3 respectively. For the π + π − cross sections shown in Fig. 2, one can find that the lines of model I, IV and V with a π 2 + a π 4 close to the same value (∼ 0.23) almost coincide with each other. And the line of model II (III) with the smallest (largest) a π 2 and a π 4 have the largest (smallest) slope † . It indicates that the twist-2 result of σ 0 (π + π − ) depends on the † The Q-dependence of the cross sections σ 0 (M + M − ) can be parameterized into the form: σ 0 (M + M − ) ∝ Q −n [11,15]. In this work, we use the "slope" to indicate the exponent n. Twist-2 results of the cross sections σ0(π + π − ) and σ0(K + K − ) with the five sample models of the distribution amplitudes listed in Table. II. The points with errors are the experimental data from Refs. [12,13,15].  Table. II. The points with errors are the experimental data from Refs. [12,13,15].
Gegenbauer momentums almost in the combination of a π 2 + a π 4 . Moreover, by fixing a π 2 (a π 4 ), i.e., comparing the result of model II with that of model IV (V), the slope of the pion lines is found to become large as a π 4 (a π 2 ) decrease. The same conclusion can also be obtained by comparing the result of model III with that of model V (IV). While, for the K + K − cross sections shown in Fig. 2, one can find that the lines of model II and V (III and IV) almost coincide with each other, which implies the twist-2 result of σ 0 (K + K − ) is insensitive to the variations of a K 1 . Furthermore, the lines of model II and V (III and IV) with the smallest (largest) a K 2 have the largest (smallest) slope, i.e., the slope of the kaon lines become large as a K 2 decrease. However, with all the five distribution amplitudes, we find that the twist-2 result of σ 0 (K + K − ) is several times smaller than its experiment data [12,13,15], and the twist-2 result of σ 0 (π + π − ) is almost an order of magnitude smaller than the data [12,13,15]. Similar conclusion was already drawn in Refs. [27,28].
At twist-3 level, to assess the uncertainties related to the Gegenbauer moments, we draw the figures at specific chiral enhancement parameters. From Fig. 2 and Fig. 3, one can see that the relative uncertainties come from the Gegenbauer moments in the twist-3 results are much smaller than that in the twist-2 results. With the Gegenbauer moments chosen at their central values and the chiral enhancement parameters being in their reasonable ranges, the twist-3 cross sections are reshown in Fig. 4. Comparing Fig. 3 with Fig. 4, it is found that the main uncertainties in the twist-3 contributions come from the chiral enhancement parameters. For both the π + π − and K + K − processes, the twist-3 cross sections are found to be in good agreement with the BELLE measurements [15] in the energy region 2.4 GeV < Q < 4.1 GeV and the ALEPH measurements [13] in the energy region 2.0 GeV < Q < 6.0 GeV. Especially, when the chiral enhancement parameters are taken the upper value µ π (2 GeV) = 2.80 GeV and µ K (2 GeV) = 2.75 GeV, the agreements are remarkable. But, in the relatively low energy region The points with errors are the experimental data from TPC [12], BELLE [15] and ALEPH [13].
1.0 GeV < Q < 2.0 GeV, there still exists considerable discrepancies between our results and the experimental data [12]. For the π + π − process, the discrepancy in 1.0 GeV < Q < 2.0 GeV may be ascribed to the interference of the continuum with the resonances, such as the f 2 (1270). For the K + K − process, the neglection of the twist-4 contributions, which are suppressed by the factor m 2 K /Q 2 , may lead to the discrepancy in 1.0 GeV < Q < 2.0 GeV. In Fig. 5, we make a comparison between the twist-2 and the twist-3 results of the cross sections σ 0 (π + π − ) and σ 0 (K + K − ) with all the nonperturbative input parameters chosen at the central values. Compared with the previous twist-2 results [19] obtained in the collinear factorization, the twist-2 results in the k T factorization are enhanced several times in the energy region 2.0 GeV < Q < 7.0 GeV, but suppressed significantly in the relatively low energy region 1.0 GeV < Q < 2.0 GeV. The enhancement in 2.0 GeV < Q < 7.0 GeV may be due to the higher-power corrections from the transverse momentum k ⊥ . While the suppression in 1.0 GeV < Q < 2.0 GeV may be ascribed to the corrections from the sudakov resummation, which have been organized into the suppression factor exp[−S]. As pointed out in Ref. [24], the asymptotic scattering amplitudes for the processes γγ → π + π − , K + K − are insensitive to the Sudakov corrections. It indicates that the Sudakov factor exp[−S] affects the cross sections slightly in the relatively high energy region, but may lead to a strong suppression at the relatively low energy scale. Comparing the twist-2 results with the twist-3 ones in the k T factorization, one can conclude that the twist-3 corrections and their interference with the k ⊥ effects are very significant for both the π + π − and K + K − processes in the intermediate energy region 1.0 GeV < Q < 7.0 GeV. Especially in the relatively low energy region 1.0 GeV < Q < 2.0 GeV, the twist-3 results of σ 0 (π + π − ) and σ 0 (K + K − ) are enhanced by almost an order of magnitude as compared with the twist-2 results. Similar situations can also be found in the calculations of the pion-photon transition form factor [50] and the pion electromagnetic form factor [37,39,46]. A more detailed investigation has shown that the twist-3 corrections fall off rapidly with increasing Q and, beyond 10 GeV region for both pion and kaon, fall below the twist-2 contributions. As is known that, the twist-2 contributions are expected to dominate in the cross section at the asymptotically large momentum transfers. However, as evident from Fig. 5, the higher twist contributions should be taken into account in the intermediate energy regions  [12], ALEPH [13] and BELLE [15].
The cross section ratio σ 0 (K + K − )/σ 0 (π + π − ) is presented in Fig. 6. We find that the twist-2 results both in the collinear factorization and in the k T factorization are almost equal to (f K /f π ) 4 ≈ 2.26. Similar theoretical results have also been given in Refs. [9,11]. While the fitted result from BELLE measurements [15] is 0.89 ± 0.04 ± 0.15 for 3.0 GeV < Q < 4.1 GeV energy region, which is clearly smaller than the twist-2 results obtained in our work and others [9,11]. Including the twist-3 corrections, our result of the ratio σ 0 (K + K − )/σ 0 (π + π − ) is found to be in line with the ALEPH measurements [13] in 2.5 GeV < Q < 4.0 GeV and the BELLE measurements [15] in 3.5 GeV < Q < 4.1 GeV. In the energy region 1.0 GeV < Q < 2.5 GeV, there still exists discrepancy between our prediction and the experimental measurements. The reason may be due to the contributions from the resonances and the m 2 M /Q 2 terms, which have been neglected in our calculations, may play a significant role in this energy region. At last but not least, we exhibit the pion and kaon angular distributions, i.e., the differential cross sections σ −1 0 dσ(γγ → M + M − )/d| cos θ|, which have been normalized to the cross section σ 0 (M + M − ). Owing to the normalization, the theoretical uncertainties of the pion and kaon angular distributions from the Gegenbauer moments and the chiral enhancement parameters cancel to a large extent. Fig. 7 shows our results for Q = 3.2 − 3.3 GeV which is in the middle range of the BELLE measurements. One can find the twist-2 results in the collinear factorization and those in the k T factorization give nearly the same angular distributions. While the angular distributions with the twist-3 corrections in the k T factorization become more steady over the region | cos θ| < 0.6. Obviously, the angular distributions obtained in this work are not in line with the BELLE measurements [15]. The reason may be related to the simplification made in Eq. (38), which is illegal in some angular coverage. In Fig. 7, we also show the known BL result [9], corresponding to the blue dot-dashed line 1.227/ sin 4 θ, which is obtained by using the approximate relation in Eq. (1) with F M (Q 2 ) ≈ 0.4 GeV 2 /Q 2 . However, the CLEO measurements [58] gave a larger form factor F M (Q 2 ) ≈ 1.01 GeV 2 /Q 2 , which would put the known BL result [9] in question. As a matter of fact, the approximate relation in Eq. (1) was obtained at the twist-2 level, and it may not hold in the intermediate energy region, where the twist-3 corrections play a key role in both the scattering processes γγ → M + M − (from Fig. 5) and the timelike electromagnetic form factor F M (Q 2 ) [38,40,50,59].

IV. SUMMARY AND CONCLUSION
Our work represents the first investigation of the twist-3 corrections to the two-photon processes γγ → π + π − , K + K − in the perturbative QCD approach based on the k T factorization theorem. The transverse momentum dependence and the resummation effects are included in the perturbative QCD approach. We use the twist-2 and twist-3 light-cone wave functions incorporating transverse degrees of freedom as the nonperturbative dynamical inputs. The nonperturbative contributions from the end-point regions can be effectively suppressed by the Sudakov factor exp[−S] and the threshold resummation factor S t (x). Consequently, the perturbative QCD calculation becomes more self-consistent and applicable, especially in the few GeV region.
Within the uncertainties from the distribution amplitudes, the twist-2 results of the cross sections σ 0 (π + π − ) and σ 0 (K + K − ) are much smaller than the experimental data [12,13,15]. As shown in our numerical analysis, it is found that both the transverse momentum effects and the twist-3 corrections play a significant role in the two-photon processes γγ → π + π − , K + K − in the intermediate energy region. In the relatively high energy regions, the twist-3 results of the cross sections σ 0 (π + π − ), σ 0 (K + K − ) and their ratio σ 0 (K + K − )/σ 0 (π + π − ) are in good agreement with their corresponding experiment data [13,15]. But in the relatively low energy region, the predicted cross sections and their ratio are still in disagreement with the experimental measurement [12]. From the above analysis, we can conclude that the cross sections for γγ → π + π − , K + K − with angular region | cos θ| < 0.6 may be dominated by the perturbative QCD contributions when the collision energy Q 2 GeV. And it is noteworthy that, by analyzing the differential cross sections, the authors in Refs. [25,26] have concluded that the transition from nonperturbative to perturbative QCD in the processes γγ → π + π − , K + K − is at Q ≈ 2 GeV and θ ≈ 40 • . Obviously, these two conclusions are consistent with each other.
However, our results of the pion and kaon angular distributions are different from the BELLE measurements [15]. The reason may be partly due to the simplification made in Eq. (38), which is illegal in some angular coverage. In addition, in the relatively low energy region 1.0 GeV < Q < 2.0 GeV, two-photon processes γγ → π + π − , K + K − have been investigated by other methods, such as the QCD sum rule [25,60] and the chiral perturbation theory [61,62], which can be seen as the complementation of the perturbative QCD approach for a full understanding of these processes.