Exclusive production of a pair of high transverse momentum photons in pion-nucleon collisions for extracting generalized parton distributions

We show that exclusive production of a pair of high transverse momentum photons in pion-nucleon collisions can be systematically studied in QCD factorization approach if the photon's transverse momentum $q_T$ with respect to the colliding pion is much greater than $\Lambda_{\rm QCD}$. We demonstrate that the leading power non-perturbative contributions to the scattering amplitudes of this exclusive process are process-independent and can be systematically factorized into universal pion's distribution amplitudes (DAs) and nucleon's generalized parton distributions (GPDs), which are convoluted with corresponding infrared safe and perturbatively calculable short-distance hard parts. The correction to this factorized expression is suppressed by powers of $1/q_T$. We demonstrate quantitatively that this new type of exclusive processes is not only complementary to existing processes for extracting GPDs, but also capable of providing an enhanced sensitivity to the dependence of both DAs and GPDs on the active parton's momentum fraction $x$. We also introduce additional, but the same type of exclusive observables to enhance our capability to explore GPDs, in particular, their $x$-dependence.


Introduction
The proton and neutron, known as nucleons, are the fundamental building blocks of all atomic nuclei, and themselves are emerged as strongly interacting and relativistic bound states of quarks and gluons of Quantum Chromodynamics (QCD). Understanding the internal structure of nucleons in terms of their constituents, quarks and gluons, and their interactions has been one of the central goals of modern particle and nuclear physics. However, owing to the color confinement of QCD, it has been an unprecedented intellectual challenge to explore and quantify the structure of nucleons without being able to see quarks and gluons directly. QCD color interaction is so strong at a typical hadronic scale O(Λ QCD ) ∼ 1/R with a typical hadron radius R ∼ 1 fm that any cross section with identified hadron(s) cannot be calculated fully in QCD perturbation theory. Fortunately, with the help of asymptotic freedom of QCD by which the color interaction becomes weaker and calculable perturbatively at short distances, the QCD factorization theorem [1] has been developed to factorize the dynamics at different momentum scales to identify good cross sections (or good physical observables) whose leading nonperturbative dynamics can be organized into universal distribution functions, while other non-perturbative contributions are shown to be suppressed by inverse power of the large momentum transfer of the collision. Predictions follow when cross sections with different hard scatterings but the same nonperturbative distributions are compared. It is the QCD factorization for physical scattering processes with a large momentum transfer Q 1/R that has enabled us to probe the particle (or partonic) nature of quarks and gluons at the short-distance, and to connect them to observed hadron(s) in terms of universal distribution functions. With a set of well determined universal distribution functions to find a quark (q), antiquark (q), or gluon (g) with a momentum fraction x inside a colliding hadron of momentum p with xp ∼ Q, known as the parton distribution functions (PDFs) f i/h (x, µ 2 ) for finding a parton of type i = q,q, g inside a colliding hadron h probed at a hard scale µ ∼ Q, QCD factorization formalism has been extremely successful in interpreting high energy experimental data from all facilities around the world, covering many orders in kinematic reach in both x and Q and as large as 15 orders of magnitude in difference in the size of observed scattering cross sections, which is a great success story of QCD and the Standard Model at high energy and has given us the confidence and the tools to discover the Higgs particle in proton-proton collisions [2,3], and to search for the new physics [4].
However, the probe with a large momentum transfer Q ( 1/R) is so localized in space that it is not very sensitive to the details of confined three-dimensional (3D) internal structure of the colliding hadron, in which a confined parton should have a characteristic transverse momentum scale k T ∼ 1/R Q and an uncertainty in transverse position b T ∼ R 1/Q. Recently, new and more precise data are becoming available for twoscale observables with a hard scale Q to localize the collision to probe the partonic nature of quarks and gluons along with a soft scale to be sensitive to the dynamics taking place at O(1/R). In addition, theoretical advances over the past decades have resulted in the development of QCD factorization formalism for two types of two-scale observables, distinguished by their inclusive or exclusive nature, which enables quantitative matching between the measurements of such two-scale observables and the 3D internal partonic structure of a colliding hadron. For inclusive two-scale observables, one well-studied example is the production of a massive boson that decays into a pair of measured leptons in hadron-hadron collisions (known as the Drell-Yan process), as a function of the pair's invariant mass Q and transverse momentum Q T in the Lab frame [5]. When Q 1/R, the production is dominated by the annihilation of one active parton from one colliding hadron with another active parton from the other colliding hadron, including quark-antiquark annihilation to a vector boson (γ, W/Z) or gluon-gluon fusion to a Higgs particle. When Q Q T 1/R, the measured transverse momentum of the pair is sensitive to the transverse momenta of the two colliding partons before they annihilate into the massive boson, providing the opportunity to extract the information on the active parton's transverse motion inside the colliding hadron, which is encoded in transverse momentum dependent (TMD) PDFs (or simply, TMDs), f i/h (x, k T , µ 2 ) [6]. Like PDFs, TMDs are universal distribution functions to find a quark (or gluon) with a momentum fraction x and transverse momentum k T from a colliding hadron of momentum p with xp ∼ µ ∼ Q k T , and describe the 3D motion of this active parton, its flavor dependence and its correlation with the property of the colliding hadron, such as its spin [7][8][9][10][11]. However, the probed transverse momentum k T of the active parton in the hard collision is not the same as the intrinsic or confined transverse momentum of the same parton inside a bound hadron. When the colliding hadron is broken by the large momentum transfer of the collision, a parton shower (the collision induced partonic radiation) is developed during the collision, generating additional transverse momentum to the probed active parton, which is encoded in the QCD evolution of the TMDs and could be non-perturbative, depending on the hard scale Q and the phase space available for the shower [5,12]. With more data from current and future experiments, including lepton-hadron semi-inclusive deep inelastic scatterings, better understanding of the scale dependence of TMDs could provide us with valuable information on the confined motion of quarks and gluons inside a bound hadron [13][14][15].
Without breaking the colliding hadron, the exclusive observables could provide different aspects of the hadron's internal structure. Since any cross section with identified hadron(s) cannot be calculated fully in QCD perturbation theory, it is necessary to have a hard scale Q 1/R for good exclusive observables for studying hadron's partonic structure. One classic example of exclusive hadronic observables is the high energy elastic π-scattering from atomic electrons [16], from which the electromagnetic form factor F π (Q 2 ) of the pion could be extracted as a function of the invariant mass of the exchanged virtual photon momentum q in the collision with Q 2 ≡ −q 2 ≥ 0. But, with the size and limited range of Q 2 , the extracted form factor F π (Q 2 ) did not reveal much information on the partonic nature of the pion. On the other hand, when Q 2 1/R 2 , F π (Q 2 ) could be factorized in terms of a convolution of two pion distribution amplitudes (DAs), φ π (x, µ) with momentum fraction x for an active quark, 1 − x for the corresponding antiquark and factorization scale µ, along with a perturbatively calculable short-distance coefficient function, as seen in eq. (4.1) 1 . The contributions from the pion's partonic states beyond a pair of active quark and antiquark are expected to be suppressed by powers of 1/(QR) [17]. Various experimental efforts have been devoted to measure the pion form factors at larger momentum transfers, from which the pion DAs could be extracted [18][19][20]. However, with the localized single hard interaction from the exchanged virtual photon, the factorized pion form factor F π (Q 2 ) is not very sensitive to the detailed shape of φ π (x, µ) as a function of x, other than the integral of φ π (x, µ) over x; see the discussion following eq. (4.2). Nucleon's internal structure could be much more rich and complex than the pion structure. QCD factorization of hard exclusive processes involving nucleons, such as large angle exclusive hadronic scattering, could be worked out, but, the corresponding calculations are much more difficult [17]. On the other hand, exclusive lepton-nucleon scattering with a virtual photon of invariant mass Q 2 1/R 2 could provide various two-scale observables, such as those in figure 1, where the hard scale is Q 2 ≡ −q 2 and the soft scale is t ≡ (p−p ) 2 . When Q 2 |t|, which is equivalent to requiring the time scale of the partonic hard collision ∼ 1/Q to be much shorter than the lifetime of the exchanged partonic states ∼ 1/ |t|, these two-scale exclusive processes are dominated by the exchange of an active qq or gg pair, as shown in figure 1, and can be systematically treated in QCD factorization approach [21][22][23]. The hadronic properties of the diffracted nucleon, the bottom part of the diagrams in figure 1, could be represented by generalized parton distribution functions (or simply, GPDs), f i/h (x, ξ, t, µ), where ξ ≡ (p−p ) + /(p+p ) + . The (p−p ) + = 2ξ[(p+p ) + /2] represents a total light-cone momentum transfer between the diffracted nucleon h and the hard partonic collision, where the light-cone components are defined as v ± = (v 0 ± v z )/ √ 2 for any four vector v µ . The GPDs were introduced by D. Müller et al. in 1994 [24], and their important roles in charactering hadron's partonic structure were further established by pioneering work in [25,26] and many years' theoretical development since then, which could be summarized in the reviews [27][28][29][30] and references therein.
By Fourier transforming the transverse component of the momentum transfer (p − p ) T to position space b T in the forward limit, p + → p + (or ξ → 0), the transformed GPD as a function of b T provides a transverse spatial distribution of quarks or gluons inside a colliding hadron at different values of momentum fraction x [31], That is, measuring GPDs could provide an opportunity to study QCD tomography to obtain images of the transverse spatial densities of quarks and gluons slicing at different momentum fraction x inside a colliding hadron. Their spatial b T dependence could allow us to define an effective hadron radius in terms of its quark (or gluon) spatial distributions, r q (x) (or r g (x)), as a function of x, in contrast to its electric charge radius, allowing us to ask some interesting questions, such as should r q (x) > r g (x) or vice versa, and could r g (x) saturate if x → 0, which could reveal valuable information on how quarks and gluons are bounded inside a hadron. Although we could expect that r q (x) (or r g (x)) is small at large x and increases when x decreases, as demonstrated in explicit model calculations [32], it is the precise knowledge of GPDs as functions of the parton flavor and kinematic variables, (x, ξ, t), that is needed for us to address these kinds of interesting and fundamental questions about the hadron, in particular, the proton and neutron, the fundamental building blocks of our visible world.
However, as clearly evident from the leading order diagrams in figure 1, the scattering with the exchange of a single virtual photon in figure 1 is effectively an exclusive 2 → 2 process: γ * (q) + h(p) → X(q ) + h (p ) with a final-state particle X = γ, π, J/ψ, ..., whose momentum is uniquely fixed by the virtual photon momentum q and total momentum transfer (p − p ) from the diffracted hadron (or ξ-and t-dependence of GPDs). Any sensitivity to the dependence of GPDs on the momentum fraction x, which is proportional to the relative momentum of the active quark and antiquark in figure 1(a) and (b), or the two gluons in figure 1(c), has to come from high order contribution and scale dependence of the process. More specifically, let's consider the deeply virtual Compton scattering (DVCS), first introduced in [33], as sketched in figure 1(a). The DVCS cross section can be naturally expressed in terms of Compton form factors (CFFs), which are then factorized as convolutions of GPDs with perturbatively calculable coefficients according to QCD factorization [23,26,34]. Extracting full details of GPDs from CFFs is a challenging inverse or deconvolution problem [35]. Due to the lack of sensitivity on the x-dependence for CFFs, it was shown [36] that based on a next-to-leading order analysis and a careful study of evolution effects, the reconstruction of GPDs from DVCS measurements does not possess a unique solution. Actually, two sample GPDs with different x-dependence can both fit the same CFFs [36].
Meanwhile, new exclusive diffractive processes have been introduced to enhance our capability to extract various GPDs from experimental measurements. Instead of the leptonnucleon scattering in figure 1, it was proposed to study the diffractive photo-production of a massive photon pair: γ(q) + N (p 1 ) → γ(k 1 ) + γ(k 2 ) + N (p 2 ) with the pair's invariant mass M γγ Λ QCD [37][38][39][40]. Similarly, the diffractive photo-production of a massive photon and meson pair: [42], as well as the diffractive production of two jets with a large invariant mass [43][44][45] were also proposed. Unlike the lepton-scattering processes in figure 1, whose factorization was proved by Collins et al. [21][22][23], the challenge for these new processes has been the lack of the same level of justification for the QCD factorization.
In this paper, we study exclusive pion-nucleon diffractive production of a pair of high transverse momentum photons: π(p π ) + N (p) → γ(q 1 ) + γ(q 2 ) + N (p ), as sketched in figure 2(a), with the photon's transverse momentum with respect to the collision axis between the colliding pion and the quark-antiquark pair from the diffracted nucleon being q T = |q 1T | = |q 2T | Λ QCD . Similar to the exclusive Drell-Yan process in pion-nucleon collision in figure 2(b) [46], or the exclusive deeply virtual lepton-hadron scattering processes in figure 1, the πN scattering process of our consideration in figure 2(a) is also a 2 → 3 exclusive process with a diffractive nucleon. Instead of measuring the lepton pair from the decay of a massive virtual photon in figure 2(b), or the scattered lepton to have the deeply virtual photon in figure 1, the hard scale of this new type of exclusive two-scale processes is provided by the large transverse momentum q T , which flows between the two back-to-back photons. The soft scale of this new type of two-scale processes is provided by t = (p − p ) 2 , the invariant mass squared of momentum transfer from the diffractive nucleon, which is the same as the soft scale of those exclusive processes in figure 1 and 2 we demonstrate that this new observable can be systematically studied in terms of QCD factorization approach with the same level of justification as those in figure 1, and our factorization arguments can be generalized to the similar type of exclusive processes, including some mentioned above. We also show that this observable can be not only a good probe of the factorized GPDs, complementary to those known exclusive processes, but also capable of providing more sensitivity to the much needed x-dependence of GPDs. When the hard scale q T is sufficiently large, the diffractive scattering on the nucleon N (p) is likely dominated by an exchange of a quark-antiquark pair, as indicated in figure 2(a), pulling more physically polarized partons into the hard collision would be suppressed by powers of q −1 T /R. Depending on the momentum flow of the active quark and antiquark, there are two distinctive kinematic regions for this exclusive process: (1) both active quark and antiquark have their momenta flowing into the hard part, as indicated in figure 2(c), and (2) only one of the active partons (quark or antiquark) has its momentum entering into the hard collision while the other has its momentum flowing out the hard collision to recombine with the spectators to form the diffracted hadron N (p ), as sketched in figure 2(a). As explained in Sec. 3, the factorization proof for these two regions requires different consideration due to the characteristic difference of soft gluons in the Glauber region. Once factorized, the region (1) gets contribution from the ERBL region of GPDs, while the other is relevant to the GPDs' DGLAP region [28]. When |t| → 0, while q 2 T Λ 2 QCD , the diffractive scattering with the nucleon in figure 2(c) is kinematically similar to the Sullivan process in lepton-nucleon scattering [47] and becomes sensitive to the nucleon's pion cloud. The production of the massive photon-pair in this kinematic regime (|t| → 0) could be viewed approximately as an annihilation of a real pion and a virtual (or almost real) pion of the colliding nucleon. To help present our justification of QCD factorization for exclusive massive photon-pair production in pion-nucleon collision in figure 2(a), we first demonstrate how the exclusive scattering amplitude of a simpler exclusive process, π + (p 1 ) + π − (p 2 ) → γ(q 1 ) + γ(q 2 ) with q T Λ QCD in figure 3, can be systematically factorized into a convolution of two pion DAs along with an infrared safe and perturbatively calculable short-distance coefficient in Sec. 2. With the large transverse momentum flow from one photon to the other through the hard scattering, interfering with the relative momentum flow between the active quark and antiquark of the colliding pion(s), the q T distribution of one of the two produced photons (or the equivalent cos θ distribution of the photon with respect to the collision axis in the pair's rest frame) can be sensitive to the momentum difference between the quark and antiquark of colliding pion(s), providing the sensitivity to the shape of factorized pion DAs.
In Sec. 3, we extend our collinear factorization arguments for the single-scale exclusive process: π + (p 1 ) + π − (p 2 ) → γ(q 1 ) + γ(q 2 ) to the two-scale exclusive observable: With the nonlocal color coherence between the incoming and the outgoing (or diffracted) nucleon, the N and N , we need additional discussions and reasoning for justifying the factorization of soft gluon interactions for this two-scale observable. We argue that when q T |t|, the leading contribution to exclusive scattering amplitude of π(p π ) + N (p) → γ(q 1 ) + γ(q 2 ) + N (p ) can be factorized into the universal GPDs convoluted with a pion DA along with infrared safe and perturbatively calculable coefficients. The corrections to this factorized expression is suppressed by powers of |t|/q 2 T . We show that by extending the π + π − process to πN process, the scattering amplitude develops both real and imaginary parts, both of which contribute to the cross section, and contains contributions from both unpolarized and polarized GPDs. Consequently, this new type of two-scale exclusive processes can be sensitive to both unpolarized and polarized GPDs.
In Sec. 4, we demonstrate numerically the sensitivity of this new type of exclusive high transverse momentum observables to the functional forms of pion DAs and nucleon GPDs in terms of their x-dependence. We introduce a flexible parametrization for DAs and a simplified version of the GK model for nucleon GPDs [48][49][50] with parameters to adjust their dependence on the parton momentum fraction. With our perturbatively calculated short-distance coefficients and our models for nucleon GPDs and pion DAs, we show explicitly how sensitive this exclusive production of a pair of high-q T photons can be to the shape of nucleon GPDs and pion DAs as functions of x. We also point out that such sensitivity could be enhanced with improved high-order calculation of the short-distance coefficients so that they are more perturbatively reliable at the end points where the momentum fraction of active parton from DAs and GPDs vanishes. Finally, in Sec. 5, we present our summary and outlook on opportunities to measure this new type of exclusive process at J-PARC and other facilities. We also discuss possibilities of additional two-scale observables of this type, which have the hard scale provided by the large transverse momentum q T of two exclusively produced "back-to-back" final-state particles (or jets) with q T |t| Λ QCD . The results of the hard coefficients are presented in the Appendix.
2 Exclusive production of a pair of high transverse momentum photons in a π + π − annihilation Exclusive production of a pair of high transverse momentum photons in π + π − annihilation, as sketched in figure 3, has a single observed hard scale, q T , the transverse momentum of one of the two produced photons with respect to the π + π − collision axis. The large scale q T leads to a point-like interaction that is sensitive to the partonic structure of the pions. It is then natural to consider QCD collinear factorization approach for studying this exclusive process. We show in this section that when q T Λ QCD , the scattering amplitude of this exclusive process can be factorized in terms of two pion DAs and a perturbatively calculable hard part, with corrections suppressed by powers of 1/q T . One of the main steps in deriving the factorization is to deform the soft gluon momenta out of the Glauber region. This is straightforward for the π + π − annihilation process because there is no pinch in the Glauber region, as we will show below. When we generalize the factorization formalism to the diffractive πN process in Sec. 3, an additional kinematic region, referred to as DGLAP region for GPD, appears, for which the soft gluon momentum is partly pinched in the Glauber region, and some modification is needed to prove the factorization.
figure 3: (a) Exclusive massive photon-pair production in π + π − annihilation, and (b) sample diagram to show the existence of perturbative pinch needed to separate the dynamics at different scales.

The process and corresponding kinematics
We study the exclusive production of a pair of high transverse momentum back-to-back photons in π + π − annihilation in the center-of-mass (CM) frame of the collision, as sketched in figure 3(a), where π + moves along +ẑ direction and π − along −ẑ direction. The scattering amplitude of this exclusive process is defined as where ε λ µ (q) is the polarization vector for a photon of momentum q and polarization λ. In the CM frame of this 2 → 2 process, the energy of the colliding pion is the same as the energy of the observed photon and equal to √ s/2 with s = (p 1 + p 2 ) 2 = (q 1 + q 2 ) 2 = s γγ Λ 2 QCD . By requiring q T Λ QCD , we can safely neglect the pion mass m π in the following discussion of the leading power QCD factorization of this process in the power expansion of 1/q T .

All-order factorization of exclusive scattering amplitude
The development of factorized cross sections starts with an examination of scattering amplitudes in terms of general properties of Feynman diagrams in QCD perturbation theory. When q T ∼ √ s becomes large, the exclusive π + π − annihilation process, as sketched in figure 3(a), is associated with two distinctive scales: (1) the hard scale q T characterizing the short-distance (perturbative) hard collision to produce the massive photon pair, as shown by the middle blob of the diagram in figure 3(b), and (2) the soft scale O(m π ) ∼ Λ QCD characterizing the long-distance (non-perturbative) hadronic dynamics associated with the colliding pions. A consistent separation of QCD dynamics taking place at these two distinctive scales can lead to a factorization formalism, which is an approximation up to corrections suppressed in power of m π /q T . The validity of perturbative QCD factorization formalism requires the suppression of quantum interference between the dynamics taking place at these two different momentum scales. That is, the dominant contributions to the factorized formalism should necessarily come from the phase space where the active parton(s) linking the dynamics at two different scales are forced onto their mass shells, and are consequently long-lived compared to the time scale of the hard collision. For example, for the exclusive scattering amplitude in figure 3(b), the suppression of quantum interference between the dynamics taking place in the middle blob at O(q T ) and the blobs on its left and right associated with the colliding pions requires us to demonstrate that the active quark, antiquark and gluon(s) from the colliding pions are effectively forced to be near their mass shells. However, all internal loop momentum integrals to any scattering amplitude are defined by contours in complex momentum space, and it is only at momentum configurations where some subset of loop momenta are pinched that the contours are forced to or near mass-shell poles that correspond to long-distance behavior. These "pinch surfaces" in multidimensional momentum space can be classified according to their reduced diagrams, found by contracting off-shell lines to points, from which we then derive the factorization formalism.

Reduced diagrams and leading pinch surfaces
Reduced diagrams specify the regions in the multidimensional loop momentum space that give dominant contributions to the loop integrals. Such leading regions are more conveniently realized in cut diagram notation of inclusive cross sections, in which graphical contributions to the cross sections are represented by the scattering amplitude to the left of the final state cut and the complex conjugate amplitude to the right. In the complex conjugate graphs all roles of momentum integrals are reversed with an opposite sign of i , which are responsible for the pinched poles associated with initial-and final-state interactions. However, for the factorization of exclusive scattering amplitudes, like the one in figure 3(b), all partons are internal and virtual. Their pinched poles, if there is any, do not come from the pair of the same propagators in the amplitude and its complex conjugate amplitude, since the momentum flows through them do not have to be the same in the amplitude and its complex conjugate amplitude. For the exclusive scattering amplitudes, like the one in figure 3(b), it is the integration of the relative momentum of any two active partons that pinches their momenta to be approximately on mass-shell if the invariant mass of these two active partons from the colliding pion is much smaller than their total energy.
We illustrate this pinch of loop momenta by using the sample diagram in figure 3(b) and labeling the active quark and antiquark momenta from π + on the left as k q = K/2 + k and kq = K/2 − k, respectively. The scattering amplitude in figure 3(b) then takes the form, whereĤ andR π − represent the middle blob and right-hand-side of the diagram, respectively, the ⊗ l j and ⊗ k i indicate the convolution of parton momenta l j and k i , respectively, with i, j = 1, 2, ...,D π + represents the DA of the π + of momentum p 1 , and K and k are the total and relative momentum of the active quark-antiquark pair on the left. If the total momentum of the pair is dominated by K + , we can identify the relevant perturbative contribution from the integration of k in eq. (2.6) by examining the pole structure of its k − integration. From the denominators of eq. (2.6), we have the two poles for k − , where we neglected the quark mass and overall transverse momentum of the pair K T . These two denominators pinch the k − integral, when the total energy of the pair (or its light-cone momentum K + ) is much larger than the virtuality of the pair, so long as we are away from the region k + → ±K + /2, where the quark (or antiquark) of the pair carries all the momentum while the other carries none. We should assume that this region is strongly suppressed by the π's DA when p + 1 m π . It is then clear from eq. (2.7) that the contributions from the diagram in figure 3(b) are forced into the region of phase space where the active quark and antiquark are both close to their mass shells. The same consideration can be applied to any pair of almost parallel active partons from the nonperturbative blob either on the left or the right in figure 3(b). That is, at the amplitude level, pinches happen among each pair of collinear partons from either the π + or π − side, as long as their total energy (or light-cone plus or minus momentum) is much greater than their invariant mass, which means that those partons evolve well before they enter the short-time hard interaction. Therefore, it is possible to factorize the two non-perturbative blobs associated with π + and π − , respectively, from the short-distance hard scattering process.
The generalization of the above pinch analysis leads to the so-called Libby-Sterman analysis [51,52], by which all loop momenta can be categorized into three groups: hard, collinear and soft, which we do not repeat here. Each external particle is associated with a group of collinear lines. With the assumption that the two observed high-q T photons are produced in the same hard scattering, the relevant reduced diagrams for the exclusive π + π − → γγ scattering amplitude are illustrated in figure 4. At the pinch surfaces, there are two groups of collinear lines associated with the directions of colliding π + and π − , respectively, shown as solid lines, and a hard part for the exclusive production of two backto-back high-q T photons. We can have an arbitrary number of collinear lines attaching the collinear subgraphD 1 orD 2 to the hard subgraphĤ. In addition, we can have an figure 4: General reduced diagram for the scattering amplitude of exclusive annihilation process for π + π − → γγ. Both dashed and solid lines can represent quarks and/or gluons.
arbitrary number of soft lines attaching toD 1 ,D 2 and/orĤ, represented by the dashed lines. Important contributions to the exclusive scattering amplitude come from the neighborhood of the pinch surfaces characterized by the reduced diagram in figure 4, but, not all of them contribute to the leading power term in 1/q T expansion. The leading pinch surfaces, contributing to the leading power, can be identified and determined by performing a power-counting analysis for the neighborhood of the reduced diagram in figure 4.
We characterize these regions of momentum space by introducing dimensionless scaling variables, denoted as λ, which control the relative rates at which components of loop momenta vanish near the pinch surfaces. A leading region is the one for which a vanishing region of loop momentum space near a pinch surface produces leading power behavior in 1/q T expansion. For the loop momenta k i and l j , which attachD 1 andD 2 to theĤ, respectively, we choose k i ∼ 1, λ 2 , λ Q and l j ∼ λ 2 , 1, λ Q , (2.8) with Q ∼ O(q T ) Λ QCD as a characteristic hard scale and λ ∼ O(Λ QCD /Q). We have k 2 i ∼ O(λ 2 Q 2 ) → 0 and l 2 j ∼ O(λ 2 Q 2 ) → 0 to quantify how the loop momenta approach to the pinch surface as λ → 0. We choose the momentum of the soft loops to have the following scaling behavior, with all components vanishing at the same rate, maintaining k 2 s ∼ O(λ 2 s Q 2 ) → 0. In principle, the two scaling variables, λ and λ s , need not be the same or related. In our discussion of power-counting, we choose λ s ∼ λ 2 . Considering the sample diagram in figure 5, we have That is, for the leading power contribution, we only need to keep the components k − s and k + s for soft gluon momentum k s to enter theD 1 andD 2 , respectively. A more comprehensive discussion including power-counting for subdivergences, as some loop lines approach the mass-shell faster than others, can be found in [53]. Following the power-counting analysis arguments in [6,22], we obtain the scaling behavior for the reduced diagram in figure 4 as where ⊗ indicates both the contraction of Lorentz indices, spinor indices and convolution of loop momenta, and the power α is given by 12) where N q (A→B) , N g (A→B) and N g (A→B) represent, respectively, the number of quarks, gluons and physically polarized gluons connecting from subgraph A to B. It is clear from eqs. (2.11) and (2.12) that the leading pinch surfaces (or leading regions) to the scattering amplitude are those in figure 6 with the minimal power α = 2. Given the fact that the meson π ± has one valence quark and antiquark, we must have a pair of quark and antiquark lines out of bothD 1 andD 2 for this exclusive scattering process in figure 6. At α = 2, all the gluons linkingD 1 (andD 2 ) toĤ (and S) are longitudinally polarized. Although the pinch surfaces in figure 6b are expected to provide the leading contributions from the perturbative power-counting analysis, reasonable arguments would make these contributions power suppressed. One simple argument is to recognize that with the quark lines from the soft factor S, these contributions are likely proportional to the end point of the pion wave function, where one of the two valence quarks carries almost no momentum. Since the pion wave function is expected to vanish at this point, we could conclude that figure 6b does not contribute at the leading power, but, might impact factorization at higher powers. Another possible argument for the contribution from figure 6b to be power suppressed could be achieved by studying the situation in which the soft loop momentum is scaled with λ s ∼ λ [22]. figure 6: Two possible leading regions for the exclusive π + π − annihilation to two photons. An arbitrary number of longitudinally polarized gluons can connect the collinear subgraphD 1 orD 2 toĤ or S.

Approximations
With the leading region identified as figure 6a, we introduce some controllable approximations to pick up the leading power contributions from the Feynman diagrams to prepare ourselves for performing the factorization of the collinear and soft gluons in next two subsections, respectively.
We first introduce two sets of auxiliary vectors to help extract leading contributions from the collinear and soft regions, respectively, where the non-light-like vectors n 1 and n 2 are introduced to regulate rapidity divergence with finite parameters y 1 and y 2 to keep them slightly off light cone. To avoid confusion of notations, in this subsection, we do not use the n andn as in eq. (2.4).
The leading reduced diagram in figure 6a can be formally expressed and approximated at the leading power as, where k q and kq (l q and lq) are the active quark and antiquark momenta from the collinear partD 1 (D 2 ) to the hard partĤ, respectively, {k} = k 1 , k 2 , ...
and similarly, In deriving eq. (2.15), we made the following approximations for all parton momenta to pick up leading power contribution, In eqs. (2.14) and (2.15), corresponding convolution of loop momenta (or momentum components) are suppressed. In deriving the first three lines of eq. (2.15), we used to pick up gluons' Lorentz components that give the leading power contribution in the Feynman gauge, where we suppressed the color indices. In eq. (2.19), the iε prescription is chosen such that the poles of k i · w 2 = k + i and l j · w 1 = l − j introduced by the inserted factors do not affect the deformation of soft Glauber gluon momentum discussed later. Even though we are considering the collinear gluons now, the same momenta can also reach soft region and especially the Glauber region, which are treated coherently, and the approximators must be applied on the whole momentum integration regions in order for the use of subtraction formalism for the overlapping regions [6]. In deriving the last line of eq. (2.15), we used where the color indices are again suppressed and the sign of iε as well as the space-like n 1,2 vectors in eq. (2.20) are chosen to be compatible with the contour deformation of "+" and "−" components of soft momenta when they are in the Glauber region as discussed later.
With the large {k + } inD 1 and {l − } inD 2 , we only need to keep {k − s } inD 1 and {l + s } in D 2 , respectively, in eq. (2.15) for ensuring the leading power contributions, as indicated in eq. (2.10) 2 . This is justified for the canonical scaling in eq. (2.9), but may not be valid for the soft momenta k s in the Glauber region, where we have soft gluons with figure 5 as an example, where the propagator (k i + k Glauber and (k Glauber sT ) 2 terms. As part of the soft region that also gives leading-power contribution, the Glauber region endangers factorization since it forbids the approximations made in eq. (2.15) or (2.20) which are key to the use of Ward identities (to be discussed in the next two subsections) to factorize soft gluons out of the collinear subgraphsD 1 and D 2 .
Fortunately, in the Glauber region, we can approximate the propagator of the soft gluon of momentum k s as 1/(k 2 s + iε) → 1/(−k 2 sT + iε) to be independent of k + s and k − s . Then with neglect of k + s inD 1 and k − s inD 2 , the only poles of k + s (k − s ) come from the propagators inD 2 (D 1 ), which all lie on one side of integration contour in the complex plane because all the collinear parton lines inD 2 (D 1 ) move into the hard part with positive minus (plus) momenta, as a special feature of π + π − → γγ annihilation process. The integrations of k + s and k − s are thus not pinched in the Glauber region, so that we can get out of it by deforming the contours of integration over k + s , k − s . Specifically, for a soft gluon of momentum k s in the Glauber region to enterD 1 , we deform the k − s contour as k − s → k − s + iO(λQ), and for a soft gluon of momentum l s in the Glauber region to enterD 2 , we deform the l + s contour as l + s → l + s + iO(λQ). This deformation keeps all the components, k + s , k − s and k sT (or l + s , l − s and l sT ), effectively in the same order ∼ O(λQ), allowing us to keep only {k − s } inD 1 , and {l + s } inD 2 . Unfortunately, for the meson-baryon case, πN → γγN to be discussed in the next section, the soft gluon momentum component k − s can be trapped in the Glauber region if N is moving in the "+" direction, and additional discussion is needed for treating the Glauber region.
For extracting the leading power contribution from the spinor of active quark enterinĝ H fromD 1 or leavingĤ intoD 2 , we insert the following spinor projector [6], For a quark line enteringĤ fromD 2 or leavingĤ intoD 1 , we have corresponding spinor projector, These projectors effectively project out the largest components of the spinor indices of active quarks and antiquarks, which give the leading power contributions to the exclusive scattering amplitude in the 1/q T expansion. Among all approximations listed above, the biggest error comes from neglecting the transverse components of active parton momenta entering into H, which leads to an error of order Λ QCD /q T .
The approximate expression in eq. (2.15), with the spin projections in eqs. (2.22) and (2.23) applied to active quark and antiquark lines, represents the leading power contribution to the exclusive ππ → γγ scattering amplitude in 1/q T expansion. In next two subsections, we demonstrate that this leading power contribution can be factorized into a convolution of two universal pion distribution amplitudes with a perturbatively calculable short-distance hard part that produces the two observed high transverse momentum photons.

Ward identity and factorization of collinear gluons
With the leading power contribution to the scattering amplitude given in eq. (2.15), we can use Ward identity to factorize all collinear and longitudinally polarized gluons from the hard partĤ.  From the first line in eq. (2.15), all Lorentz indices of attached gluon lines are effectively contracted by corresponding gluon momenta,Ĥ {µ},{ν} (k q ,kq, {k};l q ,lq, {l}; q 1 , q 2 ){k µ }{l ν }, which enables the use of Ward identity. We will first consider the situation with one longitudinally polarized collinear gluon of momentum k fromD 1 toĤ, as shown in figure 7.
As an identity, summing over all the possible attachments of a longitudinally polarized gluon to theĤ is equivalent to attaching it to the active quark and antiquark lines of thê H with a minus sign, as illustrated by the first equal sign in the graphic representation in figure 7. With the spinor projectors in eqs. (2.22) and (2.23) for the active quark and antiquark lines linking theĤ andD 1 and the use of the graphic Ward identity, we can move the attachment of a longitudinally polarized gluon to an active quark (or antiquark) line of theĤ to a gauge link of the same active quark (or antiquark) of theD 1 along the direction w 2 , as illustrated by the second equal sign in the graphic representation in figure 7, multiplied by the sameĤ without the gluon attachment while its active quark (or antiquark) momentum is adjusted, where g s is the strong coupling constant and t a is the generator for the fundamental representation of SU (3) color. In order to formally factor out theĤ without attachment from D 1 in eq. (2.24b), we shifted the active quark and antiquark momenta inD 1 accordingly. And in the second line of eq. (2.24b), we also reversed the gluon momentum k such that it flows along the same direction of the gauge link. Eq. (2.24) and its graphic representation in figure 7 clearly indicate that the attachment of a longitudinally polarized gluon of momentumk fromD 1 toĤ can be effectively detached fromĤ and connected to the gauge links of active quark of momentum k q and antiquark of momentum kq along the direction of w 2 with its momentum effectively flowing through the active quark (or antiquark) line into theĤ. Similarly, by applying the Ward identity to the attachment of a longitudinally polarized gluon of momentuml fromD 2 toĤ, we can effectively detach this gluon from theĤ with its momentum effectively flowing through the active quark (or antiquark) line fromD 2 , as sketched in figure 8, where the hooks on the external quark lines ofĤ indicate the amputation with the spinor projectors in eqs. (2.22) and (2.23). The attachment of multiple collinear and longitudinally polarized gluons of momentâ k µ i i with i = 1, 2, ..., I fromD 1 (l ν j j with j = 1, 2, ..., J fromD 2 ) toĤ can be detached in the same way, by summing over all possible attachments and using the Ward identity multiple times. The corresponding factor from detaching such gluons fromD 1 (D 2 ) toĤ can be combined with the factor in front ofD 1 (D 2 ) in the second (third) line of eq. (2.15) to make up the eikonal vertices and eikonal propagators that match the expansion of a product of two gauge links in the fundamental representation along the direction of w 2 (w 1 ), one from the active quark of momentum k q (l q ) and the other to the active antiquark of momentum kq (lq), to the order with a total of I (J) gluons [54]. And then by summing over all possible numbers of attachments of collinear and longitudinally polarized gluons with I, J = 1, 2, ..., ∞, we are able to detach all of them from theĤ and attach them to two gauge links along the direction of w 2 (w 1 ), or the Wilson lines in momentum space, one from the active quark of momentum k q (l q ) and the other to the active antiquark of momentum kq (lq). That is, we are able to factorize all attachments of collinear and longitudinally polarized gluons fromD 1 (orD 2 ) to theĤ into the well-defined gauge links that become a part of theD 1 (orD 2 ), as shown in figure 8, where the red thin lines indicate the color flow between collinear subgraphs,D 1 andD 2 , and the hard subgraph,Ĥ.
As pointed out above, the Wilson lines in figure 8 are in the fundamental 3 or3 color representation, indicated by the arrows on the Wilson line which denote the color flow. With the signs of i necessitated by the deformation out of the Glauber region, we have the Wilson line as where t a is again the generator for the fundamental representation of SU(3) color and will be suppressed in the rest of this paper, and the indices i, j are color indices in fundamental representation. This Wilson line in eq. (2.25) points from x to −∞ along the direction −n µ (n 0 > 0), and is a past-pointing Wilson line, like those in parton distribution functions from factorized Drell-Yan process, which is consistent with the picture that all the collinear parton lines from colliding pions come from the past to the hard collision,Ĥ, to produce a pair of high transverse momentum photons.

Ward identity and cancellation of the soft factor
It was demonstrated in the last subsection that the collinear factorsD 1 andD 2 can be detached from the hard partĤ. But, they are still connected by soft gluons from the soft factor S, which can communicate the colors between them. In this subsection, we use the approximations in eq. (2.20), which lead to the second and third lines of eq. (2.15), and the Ward identity to decouple the soft gluon attachments betweenD 1 andD 2 to achieve the factorization that we hope to derive. As discussed in the subsection 2.2.2, we only need to consider the "−" component of the soft momentum k s flowing from S toD 1 , and "+" component of the soft momentum l s flowing from S toD 2 for the leading power contribution to the amplitude. Similar to the collinear gluons, we can apply the Ward identity in the second and third lines of eq. (2.15), respectively, trying to detach all the soft gluons from their attachments toD 1 andD 2 . However, with the Wilson lines from detaching collinear longitudinal gluons from theĤ, the collinear subgraphs,D 1 andD 2 , are more complicated. Fortunately, as in eq. (2.18c), the soft momentum k s flowing intoD 1 is approximated byk s ∝ w 2 and since the Wilson line onD 1 has the vertices proportional to w 2 , the attachment of soft gluon of momentumk s to the gauge links ofD 1 vanishes. Therefore, we are allowed to sum over all possible attachments of the soft gluons toD 1 , including to the gauge links. Consequently, the use of Ward identity allows to detach all the soft gluons that are attached toD 1 and group them into two gauge links along the direction n 1 in the same way as we detach collinear gluons from theĤ. Similarly, sincel s ∝ w 1 , we can apply the Ward identify toD 2 {ν},{σ} (l q , lq, {l}; {l s }){lσ s } to detach all the soft gluons that are attached toD 2 and group them into two gauge links along the direction n 2 on the side ofD 2 , as shown in figure 9.
With all collinear and longitudinally polarized gluons detached from the hard partĤ and their impact represented by the Wilson lines toD 1 andD 2 , as shown in figure 8, and all soft gluons detached from theD 1 andD 2 and included into gauge links to the S, as shown in figure 9, we can express the exclusive scattering amplitude for π + π − → γγ as where i, j, m, n, etc. are color indices as labeled in figure 9, the sum of repeated color indices is understood, and z 1 ≡ k + 1 /p + 1 and z 2 ≡ k − 2 /p − 2 are momentum fractions. In eq. (2.26), the soft factor is given by In deriving the collinear factorsD 1 andD 2 in eq. (2.26), we took into account the fact that at the leading power, onlyk 1 andk 2 flow into the hard partĤ µν . For a generic pion distribution amplitude (suppressing the Wilson lines), we apply the following identify, toD 1 (k 1 , p 1 ) and similar identity toD 2 (k 2 , p 2 ), and obtain the collinear factors in eq. (2.26), where T represents the time-ordering and u and d are up and down quark fields, respectively. In eq. (2.26), the spinor projectors P A and P B are given in eqs. (2.22) and (2.23), respectively, and the superscripts, µ and ν, are Lorentz indices of the two produced photons. The spinor indices in eq. (2.26), convoluting betweenĤ andD's, are suppressed, and their factorization will be discussed in the next subsection. The colliding hadrons, π + and π − in our case, are color neutral. With all the soft gluons factored out of them, the collinear factors must be in a color singlet state, (2)ij . This color contraction connects the two Wilson lines in each collinear factor to givê where the sum of repeated color indices is understood, while the spinor indices are not summed over. Φ(0, ξ − ; w 2 ) and Φ(0, ζ + ; w 1 ) are the Wilson lines in the fundamental representation, joining the u and d quark fields to make the DAs gauge invariant, which is a result of factorization. Substituting eq. (2.32) into eq. (2.26), we have where the repeated color indices, i and m are summed. The soft factor now becomes That is, the soft factor S is in fact an identity matrix in the color space, and the exclusive scattering amplitude in eq. (2.26) is effectively factorized in color space, where the repeated color indices, i and m are summed, and averaged with the factor 1/N 2 c , and the "Tr" indicates the trace over all spinor indices betweenD 1 ,D 2 , andĤ. The hard partĤ that produces the pair of high-q T photons is given by the collision of two colorsinglet, collinear, and on-shell quark-antiquark pairs, one from each colliding hadron (π + or π − in this case).
The cancellation of long-range soft gluon interactions between the colliding hadrons is essential to the factorization. It means that long-distance connections between the two collinear subgraphs are canceled, so that the evolution of each collinear part is independent. Therefore, the collinear functions can be universal. In our case, the soft gluon cancellation happens because the active parton lines entering the hard interactions are collinear and color-neutral pairs, so that the soft gluons only see a color-neutral object from each colliding hadron, and thus there are no color correlations between the two collinear systems. This is the feature for exclusive processes, which is also seen in the factorization of twoquarkonium exclusive production in e + e − annihilation [55]. But, this is different from the soft gluon cancellations for inclusive processes, for example in [56], where it is the unitarity (inclusiveness of the final states) that guarantees the soft cancellation.
We also stress that the above steps of factorizing collinear gluons and soft cancellation should be viewed as for a given order of perturbative diagram expansion with a given number of soft gluon lines andD 1 -andD 2 -collinear lines. Summing over different attachments of gluon lines in the same kinematic region amounts to summing over different diagrams with the same region decomposition, along with subtractions of smaller regions to avoid double counting. Since such subtraction does not affect the used of Ward identity, after factorizing the whole amplitude intoD 1 ,D 2 andĤ, each factor should be regarded as subtracted ones. Due to the cancellation of soft gluons, the subtracted collinear fac-torsD 1 andD 2 are the same as the unsubtracted ones in eqs. (2.33) and (2.34). And the subtracted hard factorĤ can be derived perturbatively order-by-order by using the factorization formula.

Factorization formula
After the cancellation of soft gluons, the leading power contributions to the exclusive scattering amplitude of π + π − → γγ can be factorized into the structure shown in figure 10, figure 10: The factorized form for the exclusive ππ annihilation process.
while the spinor indices fromD 1 ,D 2 andĤ are still convoluted, as indicated in eq. (2.37), which need to be disentangled.
The factorD 1 is sandwiched between P A and P B , which indicates that only the term inD 1 proportional γ − , γ 5 γ − or γ 5 γ − γ i with i = 1, 2 survives. Since π + has negative parity and zero spin, only γ 5 γ − contributes. Similarly,D 2 only has its γ 5 γ + term that contributes. The result is where the indices α, β, ρ, σ here are the spinor indices, and the distribution amplitudes for π ± are given by Charge conjugation and isospin symmetry imply D π + (z) = −D π − (z), following the convention of taking (π + , π 0 , π − ) state as an isospin triplet. This particular choice does not affect our prediction for the cross section, since it is proportional to the squared amplitude. With the separation of spinor indices, we have our final factorized expression for the exclusive π + π − annihilation amplitude, where the short-distance hard coefficient function C µν is given by where a sum over repeated indices is understood, which is effectively the trace of spinor indices. The correction to the factorized formula in eq. (2.41) is suppressed by powers of m π /q T . With the renormalization group improvement from the fact that the exclusive scattering amplitude for π + (p 1 ) + π − (p 2 ) → γ(q 1 ) + γ(q 2 ) should not depend on the specific hard scale (or, the factorization scale) at which we perform our factorization steps. And with the choice of this factorization scale, the pion distribution amplitudes and the shortdistance coefficient function in eq. (2.41) become dependent on the factorization scale µ.

Gauge invariant tensor structures for the hard coefficient
The short-distance hard coefficient, C µν (z 1 , z 2 ; p + 1 , p − 2 , q T ) in eq. (2.42), is a function of the external momenta, and its tensor structure is constrained by the symmetry of the underlying theory. The most important constraint comes from electromagnetic gauge invariance or current conservation for the two external photons which requires C µν to be expressed in terms of independent and gauge invariant tensor structures.
Because of the explicit light-cone projection γ ± in eq. (2.42), we express all external momenta, p 1 , p 2 , q 1 , q 2 , in light-cone coordinates, as we have done for q µ 1,2 in eq. (2.5). We choose three independent vectors, p µ 1 , p µ 2 and q µ T . Using q 1 · q T = −q µ T q ν T g µν = q 2 T and q 2 · q T = −q 2 T , we can write down all the independent parity-even (P-even) current-conserving tensor structures, where we defined so that with i = 1, 2. Similarly, we have parity-odd (P-odd) current-conserving tensor structures where we used the abbreviation ε µν ⊥ = ε +−µν , with the convention ε 12 ⊥ = −ε 21 ⊥ = 1. One might consider another P-odd tensor structure q 2 T ε µν ⊥ + q µ T ε ρν ⊥ q T ρ + q ν T ε µρ ⊥ q T ρ , which satisfies the current conservation in eq. (2.43). But, this tensor itself vanishes for any components of µ and ν. The next constraint is from parity conservation. If we have n γ 5 's in a given diagram, parity conservation requires corresponding scattering amplitude to satisfy which holds for each individual diagram. For π + π − scattering, n = 2 and parity conservation excludes the tensor structures in eq. (2.48). In next section, we will generalize the pion-pion process to pion-nucleon scattering, for which we will have both P-even and P-odd tensor structures. For the exclusive π + π − scattering in this section, we can express the hard coefficient in terms of a linear combination of the P-even tensors in eq. (2.45), where we introduced an overall factor −ie 2 g 2 (C F /N c )/(2s 2 ) with electric charge e, strong coupling constant g, color factor C F /N c for the leading order contribution, and collision energy squared s = 2p + 1 p − 2 to make the scalar coefficients C i = C i (z 1 , z 2 ; κ) dimensionless for i = 0, · · · , 4, with κ introduced in eq. (2.5).
Substituting eq. (2.50) in eq. (2.41), we obtain where M i is the convolution of the hard scalar coefficient C i with the normalized DA φ(z) and f π = 130 MeV being the pion decay constant.
To help simplify some long expressions in this paper, we introduce a notation We can then express our factorization formalism in eq. (2.41) as and all scalar functions M i in eq. (2.51) as with i = 0, 1, ..., 4.

The leading-order hard coefficients
At leading-order of α s , there are two types of Feynman diagrams contributing to the shortdistance hard coefficients: A) the two observed photons are radiated from the different fermion lines, which we refer to as Type-A diagrams shown in figure 11, and B) they are radiated from the same fermion line, which we refer to as Type-B diagrams shown in figure 12. With the two identical photons in the final-state, we need to consider additional diagrams that are the same as those in figure 11 and 12, but with the two photons switched. That is, we need to evaluate a total of 8 Type-A diagrams and 12 Type-B diagrams for the leading-order hard coefficients. figure 11: The Type-A diagrams, plus those with q 1 and q 2 switched. Relatively thicker gluon lines are used to indicate the large transverse momentum flow from one photon to another, in contrast to those in figure 12. figure 12: The Type-B diagrams, plus those with q 1 and q 2 switched.
In the CM frame of this exclusive scattering process, the large transverse momentum q T of one-photon should be balanced by that of the other photon. This large transverse momentum q T flows through the gluon connecting the two fermion lines for all Type-A diagrams, while it does not flow through the gluon for all Type-B diagrams. Since the relative momentum of the quark and antiquark of the colliding pion, represented by the z 1 (or z 2 ) dependence of the pion DA in eq. (2.41), flows through the gluon line of the hard scattering back to the pion, the q T -dependence of the gluon propagator of the Type-A diagrams makes the hard coefficient, and hence, the cross section dσ/dq 2 T of this exclusive process, be a sensitive probe for the z 1 (or z 2 ) dependence of pion DA. Its sensitivity depends on the relative size of contributions to the cross section from these two-types of diagrams.
For our calculation of the leading-order hard coefficients, we denote the diagrams in figure 11(a)-(d) by A1, · · · , A4, sequentially, and the ones with q 1 and q 2 switched by A1 , · · · , A4 , respectively. Their contributions to the hard coefficient C µν are denoted sequentially by C µν A1 , C µν A2 , etc. Similarly, we label the individual contribution from the Type-B diagrams in a similar way. We use C µν A and C µν B to represent total contribution from the Type-A and Type-B diagrams, respectively.
The contribution from each individual diagram in figure 11 and 12 can be evaluated by using eq. (2.42). The collinear momenta of colliding partons, as labeled in figure 11 and 12, are defined aŝ and the photon momenta are given in eq. (2.5) and also in (2.44). The external collinear quark and antiquark lines from π + on the left are contracted with γ 5 (p + 1 γ − )/2, while the collinear quark and antiquark lines from π − on the right are contracted with γ 5 (p − 2 γ + )/2. At this order, all momenta of internal propagators are completely determined. For example, for figure 11(a), we have and its contribution to C µν is given by where e u = 2/3 and e d = −1/3, s = 2p + 1 p − 2 and the observed photon momenta are defined in eq. (2.44). In eq. (2.59), the momenta of internal fermion propagators are fixed as where the dependence on the parton momentum fractions z 1,2 is factored out from the external kinematic variables. On the other hand, the exchanged gluon momentum between the two fermion lines, which is the same for the gluon propagators in all the diagrams in figure 11, is given by for which the parton fractions z 1,2 cannot be factorized out of the dependence on q T . It is this entanglement of parton momentum fractions z 1,2 and external variable q T that makes Type-A diagrams sensitive to the DA's functional form.
For the Type-B diagrams in figure 12, the momenta of the internal propagators are different, although they have the same external parton momenta. For example, we have the contribution from figure 12(a), where the momenta of internal propagators are given by Similar to eq. (2.60), all the three internal propagators, including the gluon propagator, have their dependence on parton momentum fractions z 1,2 factored out from the external kinematic variables. This is actually true for all the diagrams in figure 12. Consequently, when the hard coefficient C µν B is convoluted with DAs in eq. (2.41), the measured kinematic variables are factored out of the z 1 and z 2 integration. Therefore, for the contribution from the Type-B diagrams, changing external kinematics does not directly probe the functional form of DAs. That is, the contribution from the Type-B diagrams in figure 12 to the factorized scattering amplitude in eq. (2.41) is not directly sensitive to the functional form of DAs, but rather to their integrated values or some kind of "moments".
The contribution from one single diagram, such as that in eq. (2.59) or in eq. (2.62), is not gauge invariant, and does not have the invariant form in eq. (2.50), while the sum over all the diagrams at the same order should be gauge invariant and can be organized in the form in eq. (2.50). Actually, the sum of all Type-A diagrams and the sum of all Type-B diagrams are gauge invariant separately, and each of them can be organized into the form in eq. (2.50). We can get the contribution to the scalar coefficient C i in eq. (2.50) from each diagram by extracting the coefficient of g µν ⊥ , p µ 1 p ν 1 , p µ 2 p ν 2 , p µ 1 p ν 2 and p µ 2 p ν 1 sequentially. The terms containing q µ T or q ν T will be naturally organized such that we have the gauge-invariant form of eq. (2.50), and can be used as a cross-checking. For example, for the diagram in figure 11a, we have its contribution to each scalar coefficient C i , where we have omitted the +iε since for the simple π + π − case, those poles happen at the end points where the DA vanishes, and it does not matter to which side of the z 1 (2) integration contour the poles lie. The complete contribution to the scalar coefficients C i in eq. (2.50) from all diagrams are reorganized in compact forms and given in the Appendix, where the interplay between measured q T distribution and the z-dependence of the DA is further discussed. Charge conjugation sets some useful relations among the results of individual diagrams. Applying charge conjugation on a diagram effectively exchanges the u and d quark lines, and can be visualized by simply reversing the fermion arrows and reassigning the parton momenta such that we still have u quark carrying the momentum fraction z 1 or z 2 . At the level of calculating Feynman diagrams, this simply reverses the order of the γ matrices, which does not change the value of the Dirac trace. However, for the contraction with γ 5 γ ∓ on the π ± sides, we need to reverse γ ∓ γ 5 back to γ 5 γ ∓ , so one γ 5 would bring one extra minus sign, which leads to the following relations among the diagrams, for Type-A diagrams, and for Type-B diagrams, where we also need to exchange e u and e d . Consequently, we have the symmetry for Type-A diagrams, while for C B this symmetry is broken by the difference of e 2 u and e 2 d . These relations carry through to the scalar coefficients C i (i = 0, · · · , 4), which can also serve as a useful check; see the Appendix for some more discussion.

Exclusive differential cross section
In this paper, we focus on the exclusive production of two unpolarized back-to-back photons in π + (p 1 )π − (p 2 ) collisions, and derive the differential cross section as follows, where M 2 represents scattering amplitude squared, with final-state photon polarizations summed, where λ ε λ * µ ε λ µ = −g µµ was used. From the factorized scattering amplitude in eq. (2.51), and using we obtain the scattering amplitude square as where M i with i = 0, 1, ..., 4 can be factorized if q T Λ QCD and are given in eq. (2.56). In eq. (2.72), the terms in the bracket are dimensionless and only functions of κ = 4q 2 T /s. Therefore, from eq. (2.69) with the azimuthal angle of q T integrated, we derive the differential cross section in q T , with M 2 given in eq. (2.72). The first factor in eq. (2.73) is the Jacobian peak. The differential cross section could be smoother if one changes q T -distribution to cos θ-distribution with θ being the angle between the direction of one of the observed photon and the collision axis. One should note that the value of q T alone is not enough to completely specify an event, due to the ambiguity of whether q 1 is in the forward or backward direction of p 1 , corresponding to the ± solutions in eq. (2.5). However, since the two photons are identical, the cross section must be the same for the forward and backward configurations. Therefore, we can take the forward solution in eq. (2.73), without adding a 1/2 factor to account for the factor that two photons are identical. For the rest of this paper, we will stick to the forward solution of eq. (2.5). We can also defined the integrated "total" cross section for with q 2 T min Λ 2 QCD to ensure the factorization. In our numerical estimate below, we choose q 2 T min = 1 GeV 2 .

Exclusive production of a pair of high transverse momentum photons in diffractive meson-baryon scatterings
Having explained the main steps in factorizing the amplitude for the exclusive photon-pair production in the π + π − annihilation, we now generalize the factorization formalism to an exclusive process involving diffractive scattering of a nucleon N of momentum p, π(p 2 ) + N (p) → γ(q 1 ) + γ(q 2 ) + N (p ), (3.1) where N can be a neutron (n) or a proton (p) and π can be π + or π − , making up various exclusive processes, such as, n π + → pγγ, p π − → nγγ, p π − → Λ 0 γγ, and those that could be measured with a pion beam at J-PARC and other facilities. The pion beam can also be replaced by a kaon beam and make up more processes. As shown in figure 2(c), the exclusive process, p π − → nγγ, could be made analogous to the π + π − collision by thinking of the p → n transition as taking a virtual π + out of the proton, carrying momentum ∆ = p − p and colliding with π − to produce two hard photons exclusively. Nevertheless, the factorization cannot be trivially adapted, because that analogy only corresponds to the ERBL region of GPD for which all the active partons from the nucleon enter into the hard part and the soft gluon momentum is not pinched in the Glauber region. Now that we are considering diffractive scattering, the presence of the spectator particles from the N → N transition implies another kinematic region where some partons enter into the hard part and then come out to recombine with the spectators to form the diffracted N . This corresponds to the DGLAP region for GPD, in which the soft gluon momentum can be partly trapped in the Glauber region. Additional argument is needed for the factorization proof, which will be given in subsection 3.2.

Kinematics
In the lab frame, the nucleon (pion) is moving along +ẑ (−ẑ) direction, carrying a large plus (minus) momentum. Two photons with large and opposite transverse momenta are produced in the final state, together with a recoiled nucleon, or a baryon in general. We focus on the region of phase space where −t = −∆ 2 (q 1 + q 2 ) 2 . That is, we require that the proton be recoiled in an approximately collinear direction and the invariant mass of the momentum transfer ∆ be much smaller than the energy of this transfer. This is the condition that allows the scattering amplitude of the exclusive process in eq. (3.1) to be factorized into a transition GPD of the nucleon.
Since ∆ carries a small transverse component and sufficiently large longitudinal component, it is convenient for our analysis to boost the lab frame into the CM frame of ∆ and p 2 where ∆ is along +ẑ direction, which is also the rest frame of q 1 and q 2 , so that the ∆ could mimic the momentum p 1 of π + in the Sec. 2. We denote this frame as photon frame S γ , distinguished from the lab frame S lab . The transformation from S lab to S γ can be done by first boosting along ∆ T such that ∆ is in parallel and head-to-head with p 2 , followed by a rotation in the ∆ T -p 2 plane to make ∆ along +ẑ direction.
In the S γ frame, we have the momentum conservation, 2) and the CM collision energy squareŝ ≡ (∆ + p 2 ) 2 = (q 1 + q 2 ) 2 Λ 2 QCD . For the leading power contribution, in analog to eq. (2.3), we can parametrize ∆ and p 2 as where in the second step (and in the following) the " " means the neglect of small quantities suppressed by powers of m 2 π /Q 2 or t/Q 2 with the hard scale Q ∼ O(q T ) √ŝ . In addition, we also implicitly take the rescaled light-like ∆ and p 2 as the momenta entering the hard process, with to keep the momentum conservation manifest, which is useful for the factorization of this process.
The skewness in the lab frame is defined as where P = (p + p )/2. We then have p + lab = (1 + ξ)P + lab , p + lab = (1 − ξ)P + lab , and which defines a unique role of the skewness, quantifying the momentum flowing into the hard process from the colliding hadron of momentum p.
The invariant mass squared of the momentum transfer, t = ∆ 2 , can be related to ξ and the transverse component of the momentum transfer ∆ T by where m p is the proton mass and we neglect the mass difference between proton and neutron. For a given small t, ∆ T is bounded to be small, and ξ is effectively constrained by Every event of the exclusive process in eq. (3.1) is specified by three momenta p , q 1 and q 2 , which are constrained by on-shell conditions and momentum conservation, leading to 9 − 4 = 5 degrees of freedom in kinematics. ∆ T and ξ are sufficient to specify the neutron momentum. The photon momenta are to be described by q T in the photon frame S γ , where they are back-to-back. That is, (∆ T , ξ, q T ) fixes all the kinematics. Our process is insensitive to azimuthal angle in either ∆ T or q T , and we will integrate out these angles, leaving only three degrees of freedom, ∆ T , ξ and q T , or equivalently, t, ξ and q T as independent variables. figure 13: The general pinch-singular surface for the process (3.1).

Factorization
We generalize the factorization formula derived in Sec. 2.2 to describe the scattering amplitude of the exclusive process in eq. (3.1). As indicated by the general pinch-singular surface in figure 13, the initial-state nucleon momentum p and slightly recoiled hadron momentum p define the direction of a collinear subgraph,D 1 , which is joined by a set of collinear parton lines to the hard subgraph, from which two photons with large transverse momenta are produced. The power counting for a pinch surface is derived in the same way as what was done in the last section. The only difference is that the dimension for theD 1 is reduced by 1 because we have an extra external final-state hadron line connected toD 1 in figure 13. Like eq. (2.11), we obtain the scaling behavior for corresponding reduced diagram as where α is the same as that in eq. (2.12). With the minimum power α = 2, we obtain the leading pinch surfaces to the scattering amplitude of exclusive process in eq. (3.1), as shown in figure 14, which are slightly modified from those in figure 6. Due to the electric charge or isospin exchange,D 1 orD 2 must be connected to other subdiagrams by at least two quark lines. By the same argument at the end of Sec. 2.2.1, the pinch surface in figure 14b is power suppressed compared to that in figure 14a.

Deformation out of Glauber region
Before we adopt the approximations listed in Sec. 2.2.2 to start our factorization arguments, we note one complication that distinguishes the diffractive meson-baryon process in eq. (3.1) from the π + π − case discussed in the last section. The factorization proof of π + π − process was simplified by the fact that all the collinear parton lines go from the past to now when the hard collision takes place, without going to the future as spectators, as shown in figure 5. All the parton lines collinear to π + (π − ) have positive plus (minus) momenta, and the plus/minus momenta of the soft gluons are not trapped to be much smaller than their transverse components. We can get those soft gluons out of the Glauber region by deforming the contours of their momentum integrations, as figure 14: Two possible leading regions for the process (3.1). Arbitrary number of gluons can connect to the collinear subgraphsD 1 orD 2 from S orĤ, but they have to be longitudinally polarized. discussed in Sec. 2.2.2. However, in the πN case, or more specifically, in pπ − → nγγ case, the proton-neutron transition can have either (i) all the collinear parton lines going from the proton into the hard part, as shown in figure 15a, or (ii) some parton lines going from the proton into the hard part, but others going to the future as spectators and merging with partons coming out the hard part to form a neutron, as shown in figure 15b. The type (i) corresponds to the ERBL region of GPD, and the type (ii) is for the DGLAP region. figure 15: Difference in soft gluon interaction between ERBL region (a) and DGLAP region (b) in the elastic πN process. In (a), the k − s integration is not pinched, while in (b), the k − s integration is pinched to be in the Glauber region.
For the ERBL region, the contour deformations and approximations made to the leading regions for every possible diagram are the same as those in Sec. 2.2.2. But for the DGLAP region, the presence of proton spectator may trap the minus momenta of soft gluons at small values. For example, as shown in figure 15b, the attachment of a soft gluon of momentum k s to a spectator of the colliding proton leads to two propagators with the denominators,  figure 15b. Therefore, the argument that we used in Sec. 2.2.2 to deform the contours of plus/minus components of soft gluon momenta to get them out of the Glauber region does not work for the soft minus components in the πN case when the nucleon N is moving in the "+" direction.
Luckily, the poles for the plus components of the soft momenta are solely provided by the collinear lines from the π, which all go into the hard part with positive minus momenta. All the poles from l j + k s lie on the same half plane, and therefore, we can deform k + s as (3.11) when it lies in the Glauber region flowing into theD 2 subgraph. This is the maximal extent that we can deform k + s , which leads the soft momentum k s all the way intoD 1collinear region. That is, the soft Glauber mode is deformed to be a collinear mode, which is only possible when all the collinear lines from the π flow into the hard partĤ. Had we considered an exclusive double diffractive process: pn → pnγγ, with a pair of back-to-back high transverse momentum photons produced while the nucleons are slightly diffracted, we would have both plus and minus components of soft momenta pinched in the Glauber region, which forbids the double diffractive processes, like pn → pnγγ, pp → pp + jet + jet, etc., to be factorized into two GPDs and a hard part [57], even though there is indeed a hard scale provided by the transverse momenta of the photons or the jets.
After the deformation of Glauber gluons, we can apply all the approximations in Sec. 2.2.2. Since we will not deform k − s in DGLAP region, it does not matter what iε prescription we assign to k − s . We choose the same convention as in Sec. 2.2.2 to be compatible with ERBL region, for which we do need to deform k − s .

Soft cancellation and factorization
We first use the same arguments presented in the last section to factorize the collinear subgraphD 2 from the hard partĤ and the soft factor S. The approximation in eqs. (2.18b) and (2.19b) allows us to use Ward identity to detach all longitudinally polarized collinear gluons ofD 2 from the hard partĤ, and factorize them into Wilson lines along w 1 , as shown in figure 16. Like the π + π − → γγ case in the last section, the Wilson lines connected tô D 2 point to the past due to the choice of iε in eq. (2.19b).  figure 18: The factorized form for the process (3.1).
Next, having eqs. (2.18d) and (2.20b), we can use Ward identity to factorize soft gluons out of the collinear factorD 2 . This leaves the collinear factorD 2 uncoupled toD 1 , so that D 2 ends up being color singlet. By the same method of Sec. 2.2.4, the soft gluons coupling toD 2 cancel. The rest of the soft gluons only couple toD 1 , as in figure 17, and can be grouped intoD 1 .
We can then use eqs. (2.18a) and (2.19a), and the Ward identity to factorize all longitudinally polarized collinear gluons ofD 1 out of the hard partĤ. This step is similar to that of the π + π − case, since the soft gluon connection toD 2 has been canceled, which would have pinched the minus components of soft gluon momenta into the Glauber region. After factorizing the longitudinally polarized collinear gluons from theĤ into Wilson line, we get a color singletD 1 . Therefore, we complete the factorization arguments and have a factorized result, as shown in figure 18. The color structure of the hard part takes the same form as in eq. (2.37). But, the spinor indices are still convoluted betweenD 1 andĤ, as well as betweenD 2 andĤ, and will be dealt with in next subsection.

Factorization formula
Similar to eq. (2.37), we derived the factorized formalism for the scattering amplitude of the exclusive process in (3.1), corresponding to the factorized diagram in figure 18, where the repeated color indices, i and m are summed, and averaged with the factor 1/N 2 c , and the "Tr" indicates the trace over all spinor indices betweenD 1 ,D 2 , andĤ. In eq. (3.12), P A , P B , andD 2 (z 2 , p 2 ) are the same as those in eq. (2.37), butD 1 (z 1 , p, p ) is different, which now represents the transition GPD of the nucleon N , where α, β are spinor indices, w 2 is as in eq. (2.13), color indices have been implicitly summed, and in the second line, we shifted the position of the operator to be consistent with the convention in [28]. Now P A and P B sandwichingD 1 picks out only the term proportional γ − , γ − γ 5 or γ − γ 5 γ i . Because of helicity conservation, the transversity GPD associated with γ − γ 5 γ i does not contribute at leading power. Effectively, we have where F ud N N (z 1 , ξ, t) and F ud N N (z 1 , ξ, t) are GPDs with different chirality characterizing the amplitude for the transition of hadron N to N , where F ud N N (x, ξ, t) and F ud N N (x, ξ, t) are the GPDs defined with the convention in Ref. [28]. Note that we are using an unusual variable z 1 to label the momentum fraction of an active parton (u quark here), as indicated in figure 18, in order to have a direct analogy to the π + π − process that we studied in the last section. As clearly indicated in eqs. (3.15b) and (3.15d), the momentum fraction z 1 is closely related to the common variables of GPDs, such as x and ξ, Consequently, the range of z 1 is different from [0, 1] for the nucleon side, as opposed to z 2 for the π, and is given by The choice of z 1 parameter highlights the so-called ERBL region, which lies between −ξ < x < ξ, and is now given by 0 < z 1 < 1. In this region, a pair of quark and antiquark with positive momentum fractions enters the hard scattering. On the other hand, one of the DGLAP regions ξ < x < 1 with a quark scattering configuration corresponds to 1 < z 1 < (1 + ξ)/2ξ, while the other DGLAP region −1 < x < −ξ with an antiquark scattering configuration becomes −(1 − ξ)/2ξ < z 1 < 0. Inserting eqs. (3.14) and (2.39) into eq. (3.12) we obtain the factorized scattering amplitude for the elastic process in eq. (3.1) where C µν is the same as that in eq. (2.42) with p + 1 replaced by ∆ + , which has γ 5 attached on both proton and pion sides so is chiral even, while C µν is given by which only has one γ 5 on the pion side and is referred as chiral odd. The correction to the factorized scattering amplitude in eq. (3.18) is suppressed by an inverse power of the high transverse momentum of observed photon q T in S γ .
The hard coefficients C µν and C µν , and the factorized formalism in eq. (3.18) are manifestly invariant under a boost alongẑ. Since the transformation from S lab to S γ is only by a boost alongẑ, up to a boost and rotation characterized by ∆ T , which is neglected at leading power, the factorization formula (3.18) takes the same form in the S γ frame, and the hard coefficients C µν and C µν can be calculated in S γ , in the same way as for π + π − case.
If N = proton and N = neutron, these transition GPDs can be related to the nucleon GPDs by isospin symmetry [58] (3.20)

The leading-order hard coefficients
The leading-order diagrams are the same as those in figure 11 and 12, except that now we have two sets of hard coefficients, obtained with different spinor projectors on the nucleon side. The calculation of the chiral-even coefficients is the same as π + π − case, and the results are reorganized in a compact form in the Appendix with z 1 taking the value within [z m , z M ]. From the parity constraint (2.49), the chiral-odd coefficient C µν can be expanded into the P-odd gauge invariant tensor structures in eq. (2.48), with p 1 replaced by ∆. Similarly to eq. (2.50), we have where ∆ µ and∆ ν are defined in the same way as p µ 1 andp ν 1 in eq. (2.46), respectively. The dimensionless scalar coefficients C 1 to C 4 can be extracted from the calculated result of each diagram by using eq. (3.21), and isolating the coefficient of the term proportional to ∆ µ , p µ 2 , ∆ ν and p ν 2 sequentially. The results are collected in the Appendix. Following the discussion above eqs. (2.66) and (2.67), charge conjugation implies similar relations for the chiral-odd coefficients, but with a minus sign, i.e., for Type-A diagrams, and for Type-B diagrams. These relations carry through to each scalar coefficient C 1 , · · · , C 4 , which has been checked in the calculations. Similar to the symmetric relation in eq. (2.68), we obtain an antisymmetric relation for C µν , for Type-A diagrams, while for C B this antisymmetry is broken by the difference of e 2 u and e 2 d .

Cross section
Using eqs. (2.50), (3.21) and (3.18), we obtain the factorized scattering amplitude M µν as where with i = 0, · · · , 4 for M i or 1, · · · , 4 for M i . Like eq. (2.72), we have the full scattering amplitude squared, summing over the photon polarizations, where the average (sum) over the spins of initial-state nucleon N (final-state N ) is included in |M i | 2 and | M i | 2 . Instead of summing (or averaging) over all nucleon spins, we can introduce GPDs sensitive to the hadron spin by expressing the matrix elements of nucleon states in eq. (3.15) in terms of independent combinations of nucleon spinors and corresponding "form factors" or spin dependent GPDs, Consequently, all scattering amplitudes corresponding to independent tensor structures, M i and M i in eq. (3.27) can be expressed in terms of the spin dependent GPDs, where the superscript "[H]" means to replace the corresponding F ud N N in eq. (3.26) by H ud N N , etc. Multiplied by their complex conjugate with the spin of N (N ) averaged (summed), we have (3.31b) where the factor 1/2 for the spin average has been included. In our numerical analysis in the next section, we take |t| ≤ 0.2 GeV 2 , which constrains ξ to be ξ ≤ 0.23 by eq. (3.8). Then the terms containing M i are suppressed by a factor of about 0.1 or smaller, compared to the terms containing |M We can thus neglect them for a rough estimate. Using eq. (3.31), we can rewrite eq. (3.27) as As discussed in Sec. 3.1, we can specify an event by ∆ T , ξ and q T , with q T being the transverse momentum of the photons in the photon frame S γ . This gives dσ = 1 2s where M 2 is given in eq. (3.32) andκ = 4q 2 T /ŝ ≤ 1 is the analog of κ (defined below eq. (2.5)) for the photon system in the S γ frame. The direction of q T can be defined with respect to the N −N plane, or p-∆ T plane. But since M 2 is for unpolarized scattering, it does not depend on the azimuthal angles of q T and ∆ T , so we can integrate them out. That allows us to only use three scalars ∆ T , ξ and q T to describe the events, which by eq. (3.7) can be transformed to the three scalar variables (t, ξ, q T ), and corresponding differential cross section, dσ d|t| dξ dq 2 where B stands for the big square bracket in eq. (3.32), which is dimensionless and can be evaluated numerically once we know the pion DAs and nucleon's GPDs. In eq. (3.34), we have separated the ξ dependent factor into two parts, in which the second part, (1 + ξ)/ξ, is canceled when we integrate over q 2 T from q 2 T min toŝ/4 = ξ/(1 + ξ)(s/2).

Numerical results
In this section, we evaluate the cross sections for producing a pair of high transverse momentum photons in exclusive pion-pion and pion-nucleon scattering and test their sensitivity to the shape of DAs and GPDs in terms of active parton's momentum fraction.

End-point sensitivity and improvement from Sudakov suppression
Before we introduce our choices of DAs and GPDs to evaluate the factorized cross sections, we discuss the well-known "end-point" sensitivity associated with perturbative evaluation of factorized elastic scattering processes, and its impact on the new type of exclusive processes introduced in this paper. For the comparison, we consider the well-known perturbative calculation of pion Form Factor F π (Q 2 ), as sketched in figure 19(a), which can be extracted from elastic electronpion scattering: e( ) + π(p π ) → e( ) + π(p π ). When the momentum transfer q = − has a high virtuality, with Q 2 ≡ −q 2 Λ 2 QCD , the pion Form Factor takes the factorized form as [59], where φ is pion DA, T B (z 1 , z 2 , Q 2 ) represents the hard scattering, and the factorization scale dependence is suppressed. With the leading order diagrams in figure 19(b), the short-distance hard part is given by [59] T with color factor C F = 4/3 for SU(3) color. By substituting this lowest order hard part in eq. (4.2) into eq. (4.1), it is clear that the pion Form Factor measurement is only sensitive to the "moment" of pion DA, 1 0 dzz −1 φ(z), not the detailed shape of φ(z), even when the probing scale Q 2 varies. Although the "moment" 1 0 dzz −1 φ(z) is expected to be finite since φ(z) → 0 as z → 0, the short-distance hard part in eq. (4.2) is actually singular as z 1 (and/or z 2 ) → 0, corresponding to the situation when the virtuality of the exchanged gluon in figure 19 goes to zero and the "hard" scattering is actually not taking place at a "shortdistance". The reliability of this perturbative fixed-order calculation near the "end-point" region when z 1 (and/or z 2 ) → 0 could be improved by taking into account the "Sudakov suppression" from resumming high order Sudakov logarithmic contributions. For example, the leading order perturbatively calculated hard part in eq. (4.2) could be improved as [60] where the running coupling constant α s is evaluated at t = max( √ z 1 z 2 Q, 1/b), K 0 is the modified Bessel function of order zero and the Sudakov factor S(z 1 , z 2 , b, Q) is given in eq. (14) of Ref. [60]. In keeping the same factorized form in eq. (4.1) with the modified hard part in eq. (4.3), an evolution of pion φ(z)'s factorization scale from 1/b to the hard scale Q was neglected. With the Sudakov suppression, the perturbative hard part T B (z 1 , z 2 , Q 2 ) in eq. (4.3) is no longer singular as z 1 (and/or z 2 ) goes to zero. Like the pion Form Factor, the perturbative hard part calculated from the Type-B diagrams in figure 12 is also singular in the "end-point" region when z 1 (and/or z 2 ) → 0 or 1, as clearly evident from the behavior of the three propagators in eq. (2.64). In addition, like the hard part of pion Form Factor in eq. (4.2), the dependence on active parton momentum fractions z 1 and z 2 in eq. (2.64) is completely decoupled from the external kinematic variables, and consequently, the contribution from the Type-B diagrams to the exclusive cross section is only sensitive to the "moment" of pion DA.
On the other hand, the three propagators for the Type-A diagrams in figure 11, as shown in eqs. (2.60) and (2.61), have slightly different features. The contribution from the Type-A diagrams is less singular in the "end-point" region when z 1 or z 2 goes to zero. The dependence on active parton momentum fractions z 1 and z 2 cannot be completely decoupled from the external kinematic variables. As shown in eq. (2.61), z 1 and z 2 are entangled with externally measured photon transverse momentum q T . It is this entanglement that makes the q T -distribution of this exclusive cross section to be sensitive to the shape of the zdependence of pion DA, or GPDs in pion-baryon scattering.

Enhanced sensitivity to the shape of pion DAs
To demonstrate that the differential cross section dσ/dq 2 T for exclusive π + π − → γγ process is sensitive to both the "moment" as well as the detailed shape of pion DA, we introduce a power-form parametrization for the normalized pion DA, with α > 0 so that the "moment" 1 0 dzz −1 φ α (z) is finite. When α = 1, this normalized pion DA is effectively the same as the so-called asymptotic form of pion DA when factorization scale µ → ∞ [61]. In this subsection, we vary the power α to show how dσ/dq 2 T changes. In the following numerical calculation, we use fixed electromagnetic coupling α e = 1/137 and the one-loop running strong coupling constant α s (µ) evaluated at the scale µ = q T . For exclusive π + π − → γγ, which could be a Sullivan-type process as a part of the pπ − → nγγ diffractive scattering when the |t| is small, we choose the collision energy √ s = 3 − 6 GeV, and require q T to be greater than 1 GeV. In figure 20(a), we plot the "total" cross section defined in eq. (2.74) with q T min = 1 GeV as a function of the power α of the normalized pion DA for various collision energies. Corresponding shapes of the normalized pion DAs are shown in figure 20(b). To minimize its dependence on the collision energy, we multiplied a scaling factor s 2 to the cross section, which effectively puts all the curves with four different collision energies on top of each other. However, as shown in figure 20(a), the scaled cross section shows a very strong dependence on the value of α, which is not because the partonic hard part is a good probe of the shape of DAs. Instead, such a strong dependence on α is caused by the "end-point" sensitivity of the perturbatively calculated partonic hard part as discussed in the last subsection, and the fact, as shown in figure 20(b), that the value of pion DAs at different α have very different values near the "end-point".
Like the "Sudakov" suppression treatment for the "end-point" region of the pion Form Factor, an improvement of the "end-point" sensitivity is also needed to improve the reliability of perturbative calculation of the factorized hard parts for this new type of exclusive processes, which is beyond the scope of the current paper.
As pointed out in Sec. 4.1, the propagator of the gluon has a very different momentum structure for Type-A diagrams from those of the Type-B diagrams. The entanglement of momentum fraction z 1 , z 2 and the observed q T in the Type-A diagrams makes the q T distribution sensitive to the z-dependence of the pion DAs.
In figure 21(a), we plot the normalized q T distribution, defined as dσ/dq T divided by the total cross section σ tot ≡ σ(q T min = 1 GeV) as defined in eq. (2.74), with respect to the same normalized q T distribution evaluated with asymptotic pion DA (α = 1). Corresponding normalized pion DAs are plotted in figure 20(b). In figure 21(b), we plot the same normalized q T distribution as a function of the power α at different values of q T . The normalized q T distribution at different q T values have very different dependence on the α. Naively, from figure 21, it seems that the q T -dependence provides additional 10% sensitivity on the shape of the pion DA. Actually, the q T -dependence should have provided a much stronger sensitivity to the shape of pion DAs, if the "end-point" sensitivity of the perturbatively calculated partonic hard parts are better controlled. As pointed out in Sec. 4.1, the Type-B diagrams have a much stronger singular behavior at the "end-point" than that of Type-A diagrams. Consequently, the Type-B diagrams give a much bigger fraction of dσ/dq T from the "end-point" region of the pion DAs than what the Type-A diagrams can give, while the Type-A diagrams are more sensitive to the shape of pion DAs. If we can improve the reliability of perturbatively calculated partonic hard cross section near the "end-point" for both Type-A and Type-B diagrams, the Type-A diagrams would contribute a much bigger fraction to the differential cross section dσ/dq T , making the measurement of dσ/dq T more sensitive to the shape of the z-dependence of the pion DAs.

Enhanced sensitivity to the shape of GPDs
In this subsection, we try to demonstrate that the photon q T distribution of exclusive meson-baryon scattering process is sensitive to the functional shapes of nucleon GPD and pion DA. The dependence on pion DA was discussed in Sec. 4.2 along with the exclusive ππ annihilation process. We now focus on the sensitivity to the shape of nucleon GPD, and fix pion DA to the power-form in eq. (4.4) with α = 0.63, which is the value compatible with the Lattice QCD calculation of the second moment of DA [62]. As discussed in Sec. 3, the integration range of active momentum fraction z 1 of GPDs is extended from (0, 1) to ((ξ −1)/2ξ, (ξ +1)/2ξ), as shown in eq. (3.17), and consequently, the propagators in partonic diagrams could be on-shell leading to poles along the integration contour of z 1 . As discussed in Sec. 2.2, the reduced diagram analysis ensures that the only perturbative pinch singularity at leading power is on the lines collinear to the external hadrons, which are systematically removed from the hard part of partonic scattering and absorbed into universal long-distance DAs or GPDs. The only possible singularities of the perturbatively calculated partonic hard part could appear at the "end-point" of the integration, and need to be suppressed by the behavior of non-perturbative DAs and/or GPDs, or by improving high order perturbative calculations as discussed in Sec. 4.1. When the non-pinched pole of z 1 locates along the contour, we use the distribution identity as a practical method to deform the contour [63], where P means the principal-value integration. Our numerical integration strategy is to individually separate each pole and use Eq. (4.5) to deal with the poles on the integration contour of z 1 . In this approach, the non-pinched poles lead to imaginary parts to the scalar coefficients M i and M i of the factorized scattering amplitude, and both their real and imaginary parts contribute to the exclusive cross section through the absolute values in eq. (3.32). For our numerical analysis below, we use the kinematics of J-PARC [64] with a charged pion beam of energy around 20 GeV, as well as that of AMBER [65] with a pion beam of energy 150 GeV hitting on fixed targets. For nucleon GPD, we choose the GK parametrization [48][49][50]66], which models the GPD using double distribution, where the subscript i refers to the choice of parton flavor and nucleon GPDs H i (x, ξ, t) are defined with the convention in Ref. [28], as specified in eq. (3.15). The double distribution f i (β, α, t) is parametrized as where h i (β) is the forward PDF of flavor i, and w i is a weight function which characterizes the ξ dependence of GPD H i (x, ξ, t) and is normalized as The larger the power n i is, the less dependent H i (x, ξ, t) is on ξ. In the limit that n i → ∞, w i → δ(α), and we have which has no dependence on ξ at all. The quark double distribution is decomposed into valence and sea components, and sea quark components are taken to be the same for u sea and d sea . Since our process is only  .20)), the sea components cancel, and only valence components contribute, 3 for which we have where ε = +1 for H and −1 for H. The condition θ(β) means that H q val (x, ξ, t) = 0 only when −ξ < x ≤ 1.
In the GK model, b val = 0 for both H and H, and α val = 0.9 GeV −2 for H and 0.45 GeV −2 for H. The forward parton density h i (β) is parametrized as a "power series" of β, fitted to global-fit PDFs. It is not our purpose to use a realistic GPD, but instead we want to see how different forms of GPDs affect the q T distribution, so it is convenient to use a simple functional form for h(β), for which we choose , (4.12) which is similar to eq. (4.4) but with possibly different powers ρ and τ . The normalization factor is N = 1 for H and N = η u − η d = 1.267 for H [66]. The parameters ρ and τ are fitted to the GK model at µ = 2 GeV, and we have the best fit This gives a h(β) peaked near β = 0. We will vary the powers (ρ, τ ) around the best-fit values and compare the change of observables.

Sensitivity to GPD's x dependence
First, we examine the sensitivity of measured photon q T distribution to the x-dependence of nucleon GPDs. For simplicity we take n i → ∞ in eq. (4.8) to remove the ξ dependence : Differential cross section in eq. (3.34) as a function of photon q T for two choices of pion beam energies, along with two sets of (t, ξ) values. Different curves correspond to different (ρ, τ ) parameters for the GPD models of the chiral-even GPDs followed by that of the chiral-odd GPDs. The rise at large q T is due to the Jacobian peak of the differential cross section.
for both H and H, and have a simplified model for nucleon transition GPDs (4.14) Apart from the best-fit parameters in eq. for both H ud pn (x, ξ, t) and H ud pn (x, ξ, t). This gives a set of GPDs with their x-dependence peaked between x = 0 and 1, as shown in Figs. 22 for t = −0.1 GeV 2 . Although there is no explicit ξ dependence in eq. (4.14), the hard-part integration in (3.26) still knows about ξ since z 1 is a function of x and ξ as defined in eq. (3.16). Moreover, ξ characterizes the CM energy of the hard collision (eq. (3.6)) and thus the range of q T . Therefore, the integration of q T also differs for different ξ. As a result there will still be substantial ξ dependence of the cross section.
With the model nucleon GPDs in figure 22, we plot in figure 23 the absolute differential cross section in eq. (3.34) as a function of measured photon q T at both J-PARC and AMBER pion beam energies, along with two choices of (t, ξ) values. We have restricted q T ≥ 1 GeV to ensure that power correction to the factorization formalism is sufficiently small. The upper bound of q T depends on the collision energy and ξ. Different curves correspond to different choices of (ρ, τ ) parameters for the GPD models, which are chosen to be the same for both chiral-even and chiral-odd GPDs. The rise at large q T is due to the Jacobian peak of the differential cross section. We can avoid the Jacobian peak by plotting the differential cross sections with respect to cos θ = 1 − 4q 2 T /ŝ with θ being the angle between the observed photon and collisionẑ-axis, instead of q T , as shown in figure 24. By comparing plots on the left and right -with different ξ, and plots on the top and bottom -with increase of collision energy √ s, the q T distribution becomes more and more dominated by small q T . As √ s and ξ (or √ŝ ) increase, more phase space opens up for the production of the two back-to-back photons. As q T decreases, the virtualities of the quark propagators in the leading-order diagrams in figure 11 and 12 decrease, leading  figure 25: Ratio of normalized differential cross sections σ −1 dσ/dt dξ dq T as a function of observed photon q T evaluated with the GPD model in eq. (4.14). Different curves correspond to different parameter sets of the GPD model in eq. (4.14).
to the enhancement of differential cross sections.
To make the difference of q T shapes more manifest to better visualize the sensitivity of measured q T distribution to the x-dependence of nucleon GPDs, we plot in figure 25 the ratio of the normalized differential cross sections as a function of q T for two different collision energies, like what we plotted in figure 21. The normalized cross sections are defined by dividing the differential cross sections by σ(q T min = 1 GeV). The ratio of the normalized differential cross sections is defined by dividing by the one evaluated with the best-fit GPD model parameters in eq. (4.13) -the red curve. Taking the ratio of normalized differential cross sections effectively removes the huge variation of the absolute values of the cross sections and enhances the dependence on the parameters of GPD models, as clearly shown in figure 25. It is evident that as the peak in x-distribution of GPD model in figure 22 shifts from 0 to 1, the q T shape differs by around 10% to 20%, without even considering the possible improvement from better control of the "end-point" sensitivity as discussed in Sec. 4.1. And by comparing figure 21 and 25, we find more sensitivity to the shape of GPD than that of DA, which means the sensitivity comes more from the DGLAP region than the ERBL region. Hence, we can conclude that the shape of q T distribution has significant sensitivity to the x-dependence (or equivalently, the z 1 -dependence) of GPDs.

Sensitivity to GPD's ξ dependence
In contrast to the x-dependence of GPDs, which is proportional to the relative momentum of the active quark-antiquark pair from the diffractive nucleon, the ξ and t are direct kinematic observables once we measure the momentum of the diffracted nucleon in an event. So, in principle, getting information on ξ and t is much more direct than getting the x-dependence. However, since GPDs are collective functions of (x, ξ, t), extracting the (ξ, t) dependence of GPDs from measured (ξ, t) dependence of exclusive cross sections depends on how x-dependence is entangled with ξ-and t-dependence in GPDs, and also, in practice, how GPDs are parametrized in terms of their (x, ξ, t)-dependence.
The measured ξ-dependence of this new type of exclusive processes has three major sources: (1) ξ-dependence of GPDs, e.g., the parameter n i -dependence of the GK model in eq. (4.8); (2) ξ-dependence from the factorized scattering amplitudes, i.e., the convolution in eq. (3.26)); and (3) kinematic effect from the fact that ξ characterizes the CM energy of the hard collision when cross section is expressed in terms of (t, ξ, q 2 T ). The kinematic effect is reflected by the (1 − ξ 2 )/ξ 2 factor in eq. (3.34) and is independent of (1) and (2). In principle, it is not possible to separate the x-dependence from the ξ-dependence of GPD because of (2), i.e., the convolution of GPD and hard coefficient depends on ξ. In this subsection, we try to explore to what extent the cross section depends on how ξ is parametrized in the GPD.
To focus on the x-dependence, we set n i → ∞ in our GPD model in eq. (4.8) in our discussion in last subsection, which led to a model of GPDs that has no dependence on ξ as shown in eq. (4.14). To test the sensitivity to ξ-dependence, we choose n i = 0 and n i = 1 as two additional model GPDs. We still keep the same parametrization of h i (β) in eq. (4.12). The advantage of using small integers for n i is that we can analytically integrate out eq. (4.6) and express GPD in terms of special functions. Since our proposed process is only sensitive to the valence region, letting n i → n val , and combining eq. (4.6) with eqs. (4.7), (4.8), (4.11), and (4.12) gives us the GPD model, for n val = 1, where and is the incomplete Beta function. The parameter b v in eq. (4.7) has been set to 0. α v and N are taken unchanged from n i = ∞. In figure 26, we plot our models for chiral-even GPD H(x, ξ, t) as functions of x for three values of n val = 0, 1, ∞. Our model GPDs for n val = 0 and 1 are given in eqs. (4.16) and (4.17), respectively. For n val = ∞, the GPDs are given in eq. (4.14). We fix (ρ, τ ) to be the best-fit values (4.13), t = −0.1 GeV 2 and show GPDs for three values of ξ(= 0.1, 0.2, 0.3). GPDs with n val = 0 have the maximum ξ dependence while those with n val = ∞ have no ξ dependence, which is clearly evident from the examples of chiral-even GPD H(x, ξ, t) in figure 26.
By integrating out q T , we plot the cross section as a distribution of ξ in figure 27 for the AMBER energy E π = 150 GeV, where the kinematic factor (1 − ξ 2 )/ξ 2 has been divided out. We see that different ξ parametrizations do reflect themselves in the ξ-distribution of differential cross sections. Their relative differences are better seen by taking ratios to the one with n = ∞, as seen in figure 27(b). Comparing n = ∞ with n = 1, we see that introducing some ξ dependence to GPDs through n = 1 leads to 20% change to the ξ distribution of the cross sections. Then increasing the ξ dependence from n = 1 to n = 0 leads to a further 20% ∼ 40% change.

Sensitivity to GPD's t dependence
Same as the ξ-dependence, the t-dependence of the diffractive cross section is experimentally determined. On the other hand, the t-dependence of theoretically factorized cross section comes from (1) the t dependence of GPD and (2) kinematic effect of hard process. As shown in eq. (3.8) the value of t actually constrains the available range of the ξ. It is worth emphasizing that t does not enter the hard process directly as an immediate consequence of the leading-power factorization, which is accurate up to power corrections of |t|/Q 2 . However, the information of t is not lost, but is captured by GPDs. The Fourier transform of GPDs with respect to the transverse component of t leads to transverse spatial density distributions of quarks and gluons inside a bound hadron, which could reveal very valuable information on how quarks and gluons are distributed in an environment of a confined hadron. Comparing with the x-dependence, t-dependence is more visible in a physical process and will not be explored in more details in this work.

Discussion and Outlook
Exclusive processes provide valuable information that is different from and complementary to inclusive processes. Without breaking the hadron, exclusive diffractive processes that can be factorized into DAs and GPDs provide not only one hard scale to characterize the particle nature of quarks and gluons inside the hadron, but also a secondary soft scale t that allows us to probe into the transverse structure of the hadron to explore much needed information on spatial distributions of quarks and gluons in a bound hadron. However, the hard scale for many existing exclusive processes, such as pion Form Factor and DVCS on a nucleon, is provided by a single virtual particle, and measured exclusive cross sections are most sensitive to the total momentum transfer from the diffracted hadron, but not to the relative momentum of the quark-antiquark or two gluons from the diffractive hadron. Consequently, such exclusive processes are only sensitive to the "moments" of DAs and/or GPDs, such as 1 0 dz z −1 φ π (z), as discussed in Sec. 4. The information on such "moments" is far from enough to constrain the functional forms of DAs and GPDs, and three-dimensional spatial imaging of quarks and gluons. That is, we need to seek for more exclusive processes, like the one that we proposed in this paper, to provide better constraints on the hadron tomography.
In this paper, we introduced exclusive production of a pair of high transverse momentum photons in pion-nucleon collisions and demonstrated that this new 2 → 3 exclusive diffractive process can be systematically factorized into universal pion's distribution amplitude and nucleon's generalized parton distributions, which are convoluted with corresponding infrared safe and perturbatively calculable short-distance hard parts. The correction to the factorization of this exclusive process is suppressed by powers of Λ QCD /q T . We also showed quantitatively that this new type of exclusive processes is not only complementary to existing processes for extracting GPDs, but also capable of providing an enhanced sensitivity to the parton momentum fraction of both DAs and GPDs from the measured transverse momentum q T distribution. This new 2 → 3 exclusive process could be measured at J-PARC and AMBER. In addition, our proof of the leading-power factorization for exclusive production of a pair of high transverse momentum photons can be carried through for justifying the factorization of exclusive Drell-Yan production in πN → + − (Q)N when Q |t|, which could be measured at J-PARC and other facilities as well. We stress that the sensitivity of observed q T shape to the functional form of GPDs is because there are two back-to-back particles coming out of the hard collision and the momentum flow between these two back-to-back particles entangles with the relative momentum of the two active partons from the diffractive hadron. It is this entanglement of momenta that provides the additional sensitivity to the x-dependence of GPDs.
In addition to the exclusive production of two high transverse momentum photons in the diffractive pion-nucleon collisions, discussed in this paper, we introduce a new class of similar 2 → 3 exclusive processes for diffractive production of a back-to-back high transverse momentum pair of particles (or jets), C(p c ) and D(p d ), from a hadron h(p) to a hadron h (p ), The exclusive process in eq. (5.1) can be viewed effectively as a combination of a diffractive production of the virtual and "long-lived" partonic state(s) B * : h(p) → B * (p 2 ) + h (p ) and an exclusive production of two backto-back high transverse momentum particles (or jets) on such virtual state(s): A(p 1 ) + B * (p 2 ) → C(p c ) + D(p d ). The necessary condition for QCD factorization of the exclusive process in eq. (5.1) is that the virtuality of the intermediate state(s) B * is much smaller than its energy, i.e., the lifetime of B * is much longer than the timescale of the hard exclusive scattering to produce the two back-to-back high transverse momentum particles (or jets). It is the long lifetime of the intermediate state B * that effectively suppresses the quantum interference between the diffractive production of the B * and the hard exclusive scattering between A(p 1 ) and B * (p 2 ). This necessary condition effectively requires that the transverse momentum p c T ∼ p d T be much larger than the soft scale, We will present the sufficient condition(s) for QCD factorization of various 2 → 3 exclusive processes of the type in eq. (5.1) in a future publication. For example, from the elastic large angle pion-pion scattering, π(p 1 ) + π(p 2 ) → π(p c ) + π(p d ) we can have a new exclusive 2 → 3 diffractive production of a back-to-back pair of high transverse momentum pions: π(p 1 ) + h(p) → π(p c ) + π(p d ) + h (p ) with p c T ∼ p d T −(p − p ) 2 , if the ππ → ππ largeangle elastic scattering is dominated by a single hard scattering. Similarly, instead of the exclusive pion-baryon scattering to produce two back-to-back high transverse momentum photons, we can switch one of the final-state photon with the initial-state pion to have another 2 → 3 exclusive process: γ(p 1 ) + N (p) → γ(q 1 ) + π(p π ) + N (p ) with the back-toback high transverse momentum photon-pion pair [41,42]. Taking the advantage of photon polarization at Jefferson Lab, polarization asymmetries of an exclusive diffractive photoproduction of two high transverse momentum and back-to-back particles could provide additional channels of observables to extract various GPDs with better sensitivities on their x-dependence, which will be explored in our future publications.