Probing axial quark generalized parton distributions through exclusive photoproduction of a $\gamma\,\pi^\pm$ pair with a large invariant mass

Exclusive photoproduction of a $\gamma\,\pi^\pm$ pair in the kinematics where the pair has a large invariant mass and the final nucleon has a small transverse momentum is described in the collinear factorization framework. The scattering amplitude is calculated at leading order in $\alpha_s$ and the differential cross sections for the process are estimated in the kinematics of the JLab~12-GeV experiments. The order of magnitude of the predicted cross-sections seems sufficient for a dedicated experiment to be performed. The process turns out to be very sensitive to the axial generalized parton distribution combination $\tilde{H}_u - \tilde{H}_d$.

6. Results 13 6.1 Numerical evaluation of the scattering amplitudes and cross sections 13 6.2 Fully differential cross sections 13 6.3 Single differential cross sections 16

Introduction
Deeply virtual Compton scattering (DVCS) and deeply virtual meson production [1][2][3][4][5][6] are the two main processes under study in order to extract the generalized parton distributions (GPDs), in particular at JLab and COMPASS. The near forward photoproduction of a pair of particles with a large invariant mass is a case for a natural extension of collinear QCD factorization theorems, which allows for complementary studies of the universality of GPDs. 1 In the present paper, we extend the study of exclusive photoproduction of a γ ρ 0 pair with a large invariant mass performed by some of us [13] to the case of a γ π pair, where we limit ourselves to the production of charged pion. In both cases, two gluon intermediate state in the hard part do not contribute. The more complicated production of neutral pseudoscalar mesons is left for further studies. The process we study here is thus γ ( * ) (q) + N (p 1 ) → γ(k) + π ± (p π ) + N ′ (p 2 ) , (1.1) where (N, N ′ ) = (p, n) for the π + case and (N, N ′ ) = (n, p) for the π − case. In this process, a wide angle Compton scattering subprocess γ(qq) → γπ characterized by the large scale M γπ (the invariant mass of the final state) factorizes from generalized parton distributions (GPDs). One can relate this large scale M γπ to the large transverse momenta transmitted to the final photon and to the final meson, the pair having an overall small transverse momentum. This process is sensitive to the chiral-even GPDs due to the chiraleven character of the leading twist distribution amplitude (DA) of the pion. We believe that the experimental study of these processes should not present major difficulties to large acceptance detectors such as those developed for the 12 GeV upgrade of JLab. For the case of an outgoing pair of a charged pion and a photon, the experimental analysis should be rather easy. For the case of an outgoing pair of a neutral pion and a photon, the analysis is probably more involved since one needs to deal with a set of three photon in the final state.
Our estimated rate depends much on the magnitude of the GPDs. We will show that the expected counting rates are very sizable for a quantitative analysis, using reasonable models based on their relations to usual parton distributions.
The arguments for the factorization of fixed angle and large energy 2 → 2 processes 2 [15] allow to write the leading twist amplitude for the process γ + π → γ + π as the convolution of two mesonic distribution amplitudes and a hard scattering subprocess amplitude γ + (q +q) → γ + (q +q) with the meson states replaced by collinear quark-antiquark pairs, as illustrated in figure 1 (Left). Based on the factorization of the exclusive meson electroproduction amplitude near the forward region [16], we replace in 1 The study of such processes started in ref. [7,8] at high energy, and a similar strategy has also been advocated in ref. [9][10][11][12]. 2 The absence of any pinch singularity has been proven for cases which involves at least one photon [14]. Figure 1: Left: Factorization of the amplitude for the process γ + π → γ + π at large s and fixed angle (i.e. fixed ratio t ′ /s); Right: Replacing one DA by a GPD leads to the factorization of the amplitude for γ + N → γ + π + N ′ at large M 2 γπ .
figure 1 (Left) the lower left meson distribution amplitude by a N → N ′ GPD, and we obtain figure 1 (Right). One should note the analogy to the timelike Compton scattering process [17][18][19]: since the large lepton pair squared invariant mass Q 2 plays the role of the hard scale in a similar way as the photon-meson pair squared invariant mass for our process.
For the factorized description to apply, it is necessary to avoid the dangerous kinematical regions where a small momentum transfer is exchanged in the upper blob, namely small t ′ = (k − q) 2 or small u ′ = (p π − q) 2 , and the region where strong final state interactions between the π meson and the nucleon are dominated by resonance effects, namely where the invariant mass M 2 πN ′ = (p π + p N ′ ) 2 is not large enough. Our paper is organized as follows. In section 2, we clarify the kinematics we are interested in and set our conventions. Section 3 is devoted to the presentation of our model for DAs and GPDs. Then, in section 4, we describe the scattering amplitude of the process under study in the framework of QCD factorization, with special emphasis on the various gauge invariant classes of diagrams, which will be of importance in view of future next-to-leading studies, and the various way of fixing the gauge for the produced photon. In section 5 we explain the different steps allowing to pass from the amplitudes to the cross-sections in the most efficient way in terms of CPU time. Section 6 presents our results for the unpolarized differential cross section in the kinematics of quasi-real photon beams at JLab where S γN ∼ 6-22 GeV 2 , and we give estimates of expected rates at JLab. In appendices, we describe several technical details required by analytical and numerical aspects of our study.

Kinematics
We study the exclusive photoproduction of a meson π and a real photon on a polarized or unpolarized proton or neutron target in the kinematical regime of large invariant mass M γπ of the final photon and meson pair and small momentum transfer t = (p 2 − p 1 ) 2 between the initial and the final nucleons. Roughly speaking, in this kinematics moderate to large, and approximately opposite, transverse momenta of the final photon and meson are assumed. Our conventions are the following. We define and decompose momenta in a Sudakov basis as with p and n the light-cone vectors (2.5) The particle momenta read with M , m π the masses of the nucleon and the π meson. From these kinematical relations it follows that The total squared center-of-mass energy of the γ-N system is On the nucleon side, the squared transferred momentum is The other useful Mandelstam invariants read and The hard scale M 2 γπ is the invariant squared mass of the (γ π) system. The leading twist calculation of the hard part only involves the approximated kinematics in the generalized Bjorken limit: neglecting ∆ t in front of p t as well as hadronic masses, it amounts to For further details on kinematics, we refer to appendix A.
The typical cuts that one should apply are −t ′ , −u ′ > Λ 2 and M 2 where Λ ≫ Λ QCD and M R is a typical baryonic resonance mass. This amounts to cuts in α andᾱ at fixed M 2 γπ , which can be translated in terms of u ′ at fixed M 2 γπ and t. These conditions boil down to a safe kinematical domain (−u ′ ) min −u ′ (−u ′ ) max which we will discuss in more details in section 6. In the following, we will choose as independent kinematical variables t, u ′ , M 2 γπ . As in ref. [13], we consider here the axial gauge p µ ε µ = 0 and parametrize the polarization vector of the final photon in terms of its transverse components while the initial photon polarization is simply written as We refer to appendix B for other gauge choices, which will be relevant for future studies of next-to-leading corrections.

Non-perturbative Ingredients: DAs and GPDs
In this section, we describe the way the non-perturbative quantities which enter the scattering amplitude are parametrized.

Distribution amplitudes for the π meson
The chiral-even light-cone DA for the meson π is defined, at the leading twist 2, by the matrix element [20] π + (p π )|ū(y) with f π = 131 MeV. In the present paper, we will use the asymptotic π DA (normalized to unity)

Generalized parton distributions
In our studies, we need the p → n and n → p transition GPDs, which by isospin symmetry are identical and related to the proton GPD by the relation The chiral-even GPDs of a parton q (here q = u, d) in the proton target (λ and λ ′ are the light-cone helicities of the nucleons with the momenta p 1 and p 2 ), which are defined by [2]: We will use a parametrization of these GPDs in terms of double distributions (DDs) [21]. We refer for details to ref. [13]. In such parametrizations, GPDs are constructed from PDFs. In the present studies, we neglect any QCD evolution for these PDFs (we take a fixed factorization scale µ 2 F = 10 GeV 2 ) and we use the following models, as in ref. [13]: • For xq(x), we rely on the GRV-98 parameterization [22], as made available from the Durham database.
• For x∆q(x) , we rely on the GRSV-2000 parameterization [23], as made available from the Durham database. Two scenarios are proposed in this parameterization: the "standard", i.e. with flavor-symmetric light sea quark and antiquark distributions, and the "valence" scenario with a completely flavor-asymmetric light sea densities. We use both of them in order to evaluate the order of magnitude of the theoretical uncertainty.

The Scattering Amplitude
We now pass to the computation of the scattering amplitude of the process (2.1). When the hard scale is large enough, it is possible to study it in the framework of collinear QCD factorization, where the squared invariant mass of the (γ, π) system M 2 γπ is taken as the factorization scale. The π + meson is described as ud.
The scattering amplitude for the production of a meson π is gauge invariant, up to the well known corrections of order ∆ T √ s which have been much studied for the DVCS case [24,25]. We now concentrate on the structure of the hard part. The hard part of the diagrams is described at lowest order in α s by 20 Feynman diagrams. Half of these diagrams, denoted A and B, are drawn in figure 2. The other set (C and D diagrams) is obtained by exchanging the role of the two quarks in t−channel. This C−parity transformation corresponds to z ↔ 1 − z and x ↔ −x.

Gauge invariant decomposition of the hard amplitude
The sets of diagrams (without including charge factors) are denoted as (· · · ). We denote (AB) 123 the contribution of the sum of diagrams A 1 + A 2 + A 3 + B 1 + B 2 + B 3 , and (AB) 45 the contribution of the sum of diagrams A 4 + A 5 + B 4 + B 5 , and similarly for (CD) 12 and (CD) 345 .
They are separately QED gauge invariant. Indeed, the color factor factorizes, and the discussion reduces to a pure QED one. In the block (AB) 123 , the three bosons are connected to a single quark line in all possible ways. In the block (AB) 45 , a photon and a gluon are connected to each quark line in all possible ways. The same reasoning applies to (CD) 12 and (CD) 345 after exchanging the role of the initial and final photons.
Using the notation e q = Q q |e|, by QED gauge invariance one can write any amplitude for photon meson production as separately three gauge invariant terms, in the form where Q 1 is the charge of the quark entering the DA and Q 2 is the charge of the quark leaving the DA, in each diagram.
Considering the parity properties of the qq correlators appearing in the DA and in the GPDs, we separate the contributions for parity (+), denoted as S and parity (−), denoted as P . Only two structures occur in the hard part, namely P P (two γ 5 matrices) and SP (one γ 5 ). Now, a close inspection of the C−parity transformation which relates the two sets of 10 diagrams gives the following results. In the case of π production, for the vector contribution, the sum of diagrams reads while for the axial contribution one gets The symbol ⊗ means the integration over x (the integration over z for the pion DA is implicit and is not important here since the DA is symmetric over z ↔ 1 − z). We now denote f a GPD of the set H, E appearing in the decomposition of the vector correlator (3.4), andf a GPD of the setH,Ẽ appearing in the decomposition of the axial correlator (3.5). Let us introduce a few convenient notations. A superscript s (resp. a) refers to the symmetric (resp. antisymmetric) in x → −x part of the hard part and of the GPD, i.e.
This thus leads to and for the axial GPD contribution, i.e. P P : with Q 1 = Q u and Q 2 = Q d for a π + , and Q 1 = Q d and Q 2 = Q u for a π − . Note that this separation in QED gauge invariant blocks is somewhat simplified in the case of quarks of equal charges (π 0 or ρ 0 production), since the decompositions (4.5) and (4.6) then only involve the sum of the separately gauge invariant parts (AB) 123 and (AB) 45 (and their C-parity transforms). In the example of ρ 0 = 1 √ 2 (uū − dd) production, the amplitude is obtained by doing the simplification Q 1 = Q 2 = Q and exchanging the role of f andf . The axial GPD contribution then reads while the vector GPD contribution is (4.8) in accordance to the structure obtained in ref. [13]. One should note that contrarily to the case of ρ 0 meson production [13], which is C(−), therefore fixing a C(−) exchange in t−channel, π + production (and similarly for π − ) involves both C-parities in t−channel, which explains why both symmetrical and antisymmetrical parts of the GPDs are involved in eqs. (4.5, 4.6).

Tensor structure
For convenience, we now define the unintegrated over x and z amplitude T π through We introduce the common normalization coefficient (4.10) Note that we include the charge factors Q u and Q d inside the hard matrix element, using the decompositions obtained in eqs. (4.5, 4.6). For the P P sector, two tensor structures appear, namely Similarly, for the SP sector, the two following structures appear

Explicit computation of one diagram
As an example, we now discuss the contribution of diagram B 1 to the scattering amplitude in some details. The scattering amplitudes for π ± described by the DA (3.1) involve both the vector GPDs (3.4) and the axial GPDs (3.5). We now give the detailed expressions for T q πV [B 1 ], , for a quark with flavor q (the fact that a transition GPD is involved will be taken into account later) and for the diagram B 1 in Feynman gauge. The vector amplitude reads which includes all non trivial factors (vertices as well as quark and gluon propagators) of the hard part of diagram B 1 . The trace reads: . (4.14) Similarly one can write in the axial sector: .
For any diagram, one can now calculate its contribution to M. The integral with respect to z is trivially performed in the case of a DA expanded in the basis of Gegenbauer polynomials. The expressions for the case of an asymptotic DA, which we only consider in the present article, are given explicitly in appendix D, and expressed as linear combination of building blocks.
The integration with respect to x, for a given set of GPDs, (which can be our model described in section 3 or any other model), is then reduced to the numerical evaluation of these building block integrals.

From Amplitudes to Unpolarized Differential Cross Sections
The scattering amplitude of the process (2.1), in the factorized form, is expressed in terms of form factors H π , E π ,H π ,Ẽ π , analogous to Compton form factors in DVCS, and reads (5.1)

From amplitudes to cross sections
We isolate the tensor structures of the form factors as 3 These coefficients can be expressed in terms of the sum over diagrams of the integral of the product of their traces, of GPDs and DAs, as defined and given explicitly in appendix D. We introduce dimensionless coefficients N andÑ as follows: In order to emphasize the gauge invariant structure and to organize the numerical study, we factorize out the charge coefficients, and put an explicit index q for the flavor of the quark GPDs f q andf q . In accordance with the decompositions (4.5) and (4.6) we thus introduce 3 One should note the fact that the role of TA and TA 5 (resp. TB and TB 5 ) have been exchanged with respect to the case of γρ 0 production, see eqs. (5.13) and (5.14) of ref. [13]. This is due to the additional γ 5 structure appearing in the hard part, which can be traced through Fierz transform to the presence of a γ 5 in the matrix element (3.1).
For the specific case of our two processes, namely γπ + production on a proton and γπ − production on a neutron, taking into account the structure (3.3) of the transition GPDs structure we thus need to compute the coefficients 13) and as well asÑ Therefore, for each flavor u and d, knowing (for two given GPDs f andf , in practice H andH, see next subsection) the 12 numerical coefficients one can reconstruct the scattering amplitudes of the two processes. The expansions of these 12 coefficients in terms of 5 building block integrals are given in appendix D.2.
In this paper, we are interested in the unpolarized cross section. As a result, we will need the squared form factors after summation/average over all the polarizations (of outgoing γ and of incoming γ):

Square of M π
In the forward limit ∆ ⊥ = 0 = P ⊥ , one can show that the square of M π reads after summing over nucleon helicities: For moderately small values of ξ, this becomes: Hence we will restrict ourselves to the GPDs H,H to perform our estimates of the cross section 4 .

Cross-section
We now define the averaged amplitude squared |M π | 2 , which includes the factor 1/4 coming from the averaging over the polarizations of the initial particles. Using the expressions of the two previous subsections, and collecting all prefactors, which read 1 we have the net result, for the photoproduction of a πγ pair, Here π is either a π + or a π − , and the corresponding coefficientsÑ π + A ,Ñ π + B , N π + A 5 ,  Our central set of curves, displayed below, is obtained for S γN = 20 GeV 2 . For this value of S γN , the invariant mass M 2 γπ varies from 1.52 GeV 2 up to 9.47 GeV 2 (the crosssection vanishes at these two end points, due to the vanishing of the phase-space in −t, as shown in appendix E). We therefore vary M 2 γπ from 1.6 GeV 2 up to 9.4 GeV 2 , with a step of 0.1 GeV 2 , in order to have a full coverage of M 2 γπ for the case S γN = 20 GeV 2 . For each of these M 2 γπ values: • we calculate, for each of the above types of GPDs (in the present paper H andH), sets of u and d quarks GPDs indexed by M 2 γπ , i.e. ultimately by ξ given by The GPDs are computed as tables of 1000 values for x ranging from −1 to 1.
• we compute the building block integral I e which does not depend on −u ′ .
• we compute, for each GPD and each flavor u and d, the remaining 4 building block • this gives for each of these couples of values of (M 2 γπ , −u ′ , ), and each flavor a set of 12 coefficients listed in eqs. (5.20).

Fully differential cross sections
We now present our results for the differential cross-sections, showing in parallel the γπ + (proton target) and γπ − (neutron target) cases.     Figure 4: Left: Differential cross section for the production of a photon and a π + meson on a proton target. Right: Differential cross section for the production of a photon and a π − meson on a neutron target. Both cross-sections are at M 2 γπ = 4 GeV 2 , S γ N = 20 GeV 2 , −t = (−t) min as a function of −u ′ . In black the contributions of both vector and axial GPD, in blue the contribution of the vector GPD, and in green the contribution of the axial GPD. Solid: "valence" model, dotted: "standard" model. There is no interference between the vector and axial amplitudes.

PSfrag replacements
We first analyze the various contributions to the differential cross section in the specific kinematics: M 2 γπ = 4 GeV 2 , S γN = 20 GeV 2 , −t = (−t) min as a function of −u ′ . The dependency with respect to S γN will be discussed in section 6.4.
In figure 3, we show the relative contributions of the u− and d−quark GPDs (adding the vector and axial contributions), which interfere in a destructive way because of the flavor structure of the transition GPD. Sfrag replacements Figure 5: Left: Differential cross section for a photon and a π + meson production, for a proton target. Right: Differential cross section for a photon and a π − meson production, for a neutron target. Both are displayed as function of −u ′ , for M 2 γπ = 3, 4, 5, 6 GeV 2 (resp. in black, red, blue, green, from top to down). Solid: "valence" model, dotted: "standard" model. From this figure 3, one should note that our obtained predictions for the differential cross-sections for the production of a γπ + pair on a proton target and for the production of a γπ − pair on a neutron target are different. Indeed, contrary to a naive expectation, there is no simple relation between these two processes, since electromagnetic processes do not preserve isospin symmetry. Contrarily to the two processes γπ + → γπ + and γπ − → γπ − which are obviously C−conjugated, and thus have identical cross-sections, in our present case, the t−channel exchange is more involved. Indeed it can be interpreted as a meson exchange only in the Efremov  region −ξ < x < ξ. Technically, our processes both mix C(+) and C(−) sectors, as shown in subsection 4.1. A similar situation also occurs in the case of electroproduction of ρ + or ρ − meson, as discussed in ref. [29].
In figure 4, we show the relative contributions of the GPDs H andH involving vector and axial correlators. This demonstrates the extreme sensitivity of the differential crosssections to the axial GPDH u −H d . In the valence scenario the contributions of H u −H d and H u −H d have the same order of magnitude, while in the standard scenario, there is a clear dominance ofH u −H d . This is in contradistinction with the case of γρ 0 production [13] where the contribution of vector GPDs clearly dominates. The difference originates from the pseudo scalar nature of the pion.
We investigated the effect of changing the ansätze for the PDFs q, and thus for the GPDs H u and H d in ref. [13]. This effect was shown to be moderate, and we do not repeat this study here. Figure 5 shows the dependence on M 2 γπ . The production of the γπ pair with a large value of M 2 γπ is severely suppressed as anticipated. Note that the −u ′ range allowed by our kinematical requirements is narrower for smaller values of M 2 γπ . The two curves for each value of M 2 γπ correspond to the two parameterizations ofH(x, ξ, t), the lines corresponding to the unbroken sea scenario lying much above the other one.

Single differential cross sections
We now study the single differential cross section with respect to M 2 γπ by integrating over u ′ and t. We make a simplistic ansatz for the t−dependency of the cross-section, namely a factorized dipole form with C = 0.71 GeV 2 . The single differential cross section then reads We summarize the behavior of the domain of integration over u ′ and t when varying M 2 γπ in appendix E.

Integrated cross sections and variation with respect to S γN
For the value S γN = 20 GeV 2 , the integration over M 2 γπ of our above results within our allowed kinematical region, here 1.52 GeV 2 < M 2 γπ < 9.47 GeV 2 (see appendix E), allows to obtain the cross sections 1.2 pb < σ proton π + < 6.8 pb and 3.3 pb < σ neutron π − < 7.1 pb.  PSfrag replacements Figure 6: Left: Differential cross section dσ/dM 2 γπ + for the production of a photon and a π + meson on a proton target. Right: Differential cross section dσ/dM 2 γπ − for the production of a photon and a π − meson on a neutron target. The values of S γN vary in the set 8, 10, 12, 14, 16, 18, 20 GeV 2 . (from 8: left, brown to 20: right, blue), covering the JLab energy range. We use here the "valence"(solid) and the "standard" (dotted) scenarios.

Sfrag replacements
The variation with respect to S γN could be obtained by following the whole chain of steps described above. However, as explain in detail in ref. [13], this can be obtained almost directly from the only knowledge of the set of numerical results computed for a given value of S γN , which we take in practice as S γN = 20 GeV 2 , for any arbitrary smaller values of S γN . We summarize the idea: • we start from our set of results obtained for S γN = 20 GeV 2 , indexed by M 2 γπ and −u ′ .
• for any chosen new value ofS γN , we obtain a set of values ofM 2 γπ indexed by the set of values of M 2 γπ (which vary from 1.6 up to 9.4 GeV 2 , with a 0.1 GeV 2 step), through the relationM and for each of theseM 2 γπ a set of values of −ũ ′ , using the relation which gives the indexation of allowed values of −ũ ′ as function of known values of (−u ′ ).
As shown in ref. [13], this mapping from a given S γN to a lowerS γN provides a set of (M 2 γπ , −ũ ′ ) which exhausts the required domain. This mapping avoids the use of a very large amount of CPU time.
In figure 6 we show the differential cross section dσ/dM 2 γπ for various values of S γN covering the JLab-12 energy range. These cross sections show a maximum around M 2 γπ ≈ 2.5 GeV 2 , for most energy values. Their shapes are very similar, the only noticeable difference between the π + and the π − case being the maximum value of the differential cross-section, which is higher in the π − case. Sfrag replacements PSfrag replacements for both parameterization of the axial GPDs 5 . As for ρ 0 photoproduction, our predicted cross sections prove that the present process of photoproduction of a γπ + or γπ + pair is measurable in the typical kinematical conditions and integrated luminosity of a JLab experiment.

Counting rates
Counting rates in electron mode can be obtained using the Weizsäcker-Williams distribution. This distribution is given by [30,31] f where x is the fraction of energy lost by the incoming electron, m e is the electron mass and Q 2 max is the typical maximal value of the virtuality of the echanged photon, which we take to be 0.1 GeV 2 . We note that this distribution can be safely used based on a careful study of the scattering amplitude for the process γ * (Q 2 )N → γπ ± N ′ . This shows that in the limit Q 2 → 0, transversally polarized photons dominate and there is no appearance of any collinear singularity in this limit, which in principle could change the structure of the small Q 2 integration region [32], since in our process the quark propagators connected to the initial photon have a virtuality of the order of p 2 ⊥ ∼ M 2 γπ (this statement is valid for both ρ and π production).
Using the expression for x as a function of the incoming electron energy E e it is now easy to obtain integrated cross sections at the level of the eN process, using the relation of eq. (6.8) in figure 8.
In the case of a lepton beam, one should also consider Bethe-Heitler-type processes, in which the final real photon is emitted by the lepton beam. As discussed in ref. [13], such a Bethe-Heitler contribution is suppressed with respect to the production mechanism studied here. )d(−t)  Figure 8: Shape of the integrand of σ eN , as a function of the invariant mass of the hadronic produced state, on a proton target. Left: γπ + production on a proton target. Right: γπ − production on a neutron target. In solid-red: "valence". In dotted-blue: "standard".
The angular coverage of the final state particles is in principle a potential experimental issue. We discuss in detail the angular distribution of the outgoing photon, which might evade detection, in appendix F, taking the constraints of JLab Hall B and showing that this does not affect our predictions.
Finally, let us add a word of caution with to respect to the kinematical domain where πN ′ may enter the resonance region. A careful study of the allowed phase space shows that M 2 πN ′ (2.15) is minimal when −u ′ ∼ (−u ′ ) maxMax and M 2 γπ ∼ M 2 γπ Max . This minimal value increases with S γN . To ensure that our formalism applies, one should integrate out this domain. This is premature before precise experimental conditions are known.
We can now give our predictions for the counting rates. With an expected luminosity L = 100 pb −1 s −1 we obtain for 100 days of run: between 1.3 10 4 (valence scenario) and 8.0 10 4 γπ + pairs (standard scenario), and between 4.4 10 4 (valence scenario) and 8.9 10 4 γπ − pairs (standard scenario) in the kinematical domain discussed earlier.

Conclusion
We studied the process γN → γπ ± N ′ in the generalized Bjorken kinematics where GPD factorization is expected to hold in a collinear QCD approach. We restricted our analysis to unpolarized cross sections, which turn out to be large enough for the process to be analyzed in a quite detailed way by near-future experiments at JLab with photon beams originating from the 12 GeV electron beam.
This process is insensitive to gluon GPDs in contrast with the photoproduction of a γπ 0 pair which we leave for future studies. Our analysis has shown the dominance of the axial generalized parton distribution combinationH u −H d which is up to now not much constrained by any experimental data. Using two different reasonable ansätze based on two proposed parametrizations of polarized PDFs, we found differences by a factor of 2 to 5 in the cross-sections, see figure 6. Recent lattice studies [33] seem to favor the standard scenario which gives the larger cross section. The amplitude has very specific properties which should be very useful for future GPDs extractions programs e.g. [34].
A NLO calculation should first confirm the validity of the factorization hypothesis for this process, in the sense of absence of infrared and end-point singularities, and estimate the effects on the amplitude. Such a next to leading order computation is under study, in the spirit of ref. [35,36] in the γγ channel. Let us stress that, contrary to the DVCS (and TCS) case [19,37], the process studied here does not involve any gluonic contributions.
A similar study could be performed in the Compass experiment at CERN where S γN ∼ 200 GeV 2 and at LHC in ultraperipheral collisions [38], as discussed for the timelike Compton scattering process [39]. This also applies to future electron proton collider projects like EIC [40] and LHeC [41].
Phenomenologically, in contrast with the ρ meson, the asymptotic form of the pion DA which we have used in the present article is disputed, and quite different descriptions have been proposed. They lead to a rather universal form which was first suggested in ref. [42] and then uncovered in AdS/QCD holographic correspondence [43] as well as in dynamical chiral symmetry breaking on the light-front, see [44] and references therein. With a very good precision, it reads φ π (z) = 8 π z(1 − z) . We leave detailed studies of the impact of such types of DA for future work. As a final remark, let us stress that our study may be extended to the case of electroproduction of a photon meson pair where a moderate virtuality of the initial photon may help to access the perturbative domain with a lower value of the hard scale M γπ .

Acknowledgements
We thank Michel Fontannaz, Hervé Moutarde and Franck Sabatié for discussions. The simulations where done using the computer cluster system of CPhT. We thank the CPhT computer team for help. This

A. Some details on kinematics
In this section we give further useful expressions for kinematics.

A.1 Exact kinematics
Combining eqs. (2.11) and (2.12) one gets From eq. (2.10), we obtain so that we finally have and thus A.2 Exact kinematics for ∆ ⊥ = 0 In the case ∆ ⊥ = 0 , we now provide the exact formulas in order to get the set of parameters s, ξ, α, α π , p 2 , (−t) min as functions of M γπ , S γN , −u ′ . We refer to appendix C of ref. [13] for details. Each formula of that paper is valid here after the replacement m ρ → m π . Introducing the notationsM and Computing ξ through eq. (A.7) and then s through eq. (A.2), one can get α using The value of α π is then obtained using .
Finally, p 2 t is computed from

A.3 Approximated kinematics in the Bjorken limit
In the collinear limit, which we use for the hard part of the process,M γπ and S γN are parametrically large, and s is of the order of S γN . Neglecting ∆ 2 t , m 2 π , t and M 2 in front of s, (except in the definition of τ where we keep as usual M 2 in the denominator of eq. (A.3)), we thus have The skewedness ξ thus reads and the parameter s is given, using eq. (A.2), by

B. Electromagnetic gauge invariance
We here discuss the gauge choice for the photon polarization vectors. A first natural choice, which we also implemented in ref. [13] and use in the present article, is to consider the axial gauge p µ ε µ = 0 and parametrize the polarization vector of the final photon in terms of its transverse components while the initial photon polarization is simply written as A second choice, which will be particularly useful when computing loop corrections, is to use two different gauges for the incoming and outgoing photon to keep a symmetry between them, i.e.
One should note that as expected, ε ′µ k and ε µ k differ by a vector proportional to k, namely The expansion of ε ′µ k in Sudakov components reads so that the gauge rotation between the two Sudakov transverse components reads For ε q and ε ′ q the two transverse components are of course identical since by gauge transformation these two polarization vectors can only differ by a term proportional to q = n, in accordance to eqs. (B.2) and (B.5).

PSfrag replacements
Within the collinear factorisation approach the hard coefficient function of chiral even operators involves two types of terms: a trace over Dirac matrices without any γ 5 matrix, denoted by ǫ µ (k 1 )ǫ ν (k 2 )M µν V , and a trace over Dirac matrices with one γ 5 matrix, denoted by ǫ µ (k 1 )ǫ ν (k 2 )M µν A , with polarisation vectors ǫ's satisfying the usual orthogonality conditions ǫ(k 1 )·k 1 = 0 = ǫ(k 2 )·k 2 . Due to the momentum conservation k 1 +k 2 +P 1 +P 2 = 0 we choose as independent momenta k 1 , k 2 and P 1 . It is useful to derive the general structure of tensors M µν V and M µν A consistent with the Ward identities: The general structure of tensor M µν V

C.1.1 Decomposition of the amplitude in arbitrary gauge
In terms of the metric tensor g and independent vectors k 1 , k 2 and P 1 , due to the orthogonality conditions the general tensorial structure of M µν V can be written as where M equivalently written as One should note that only 3 of the above 4 equations are independent, so that among the 5 scalar functions M   V , which leads to the expression subject to the condition (C.8).
The gauge invariant expression for the square of the amplitude, after summing over photon polarizations, has the form or using the condition (C.8) it can be represented in the equivalent form It is useful to note that the expressions (C.11) and (C.12) correspond to results for M µν V (M V µν ) * obtained for two different gauge choices.

C.1.2 P 1 -gauge
To see that, let us first consider the light-cone gauge where k = k 1 or k = k 2 . Thus the polarization vector ǫ µ (k) having the Sudakov decomposition with respect to two light-cone vectors P 1 and k 1 can be expressed in terms of its transverse components ǫ µ ⊥ (k) satisfying the conditions ǫ µ ⊥ (k) · P 1 = 0 = ǫ µ ⊥ (k) · k 1 as In this P 1 -gauge, Using that the sum over polarisations λ equals and the expression following from the condition k 2 2 = 0, one reproduces the expression (C.11).

C.1.3 k 1 -gauge
If instead we choose the light-cone gauges ǫ(k 2 ) · k 1 = 0 (C.20) and ǫ(k 1 ) · k 2 = 0 (C. 21) which implies that in the Sudakov basis (C.13) then we reproduce directly the expression (C.12), after again using the relations (C. 18) and (C.19). The reasoning presented above is very useful for comparison of results, specially those taking into account loop radiative corrections, obtained in two different gauges.

C.2 The general structure of tensor M µν A
In a similar way the general tensorial structure of M µν A involving trace of Dirac matrices with γ 5 can be written as with A i being scalar functions constructed from the independent vectors k 1 , k 2 and P 1 . Using Schouten identity g ν µ ǫ ρστ λ + g ρ µ ǫ στ λν + g σ µ ǫ τ λνρ + g τ µ ǫ λνρσ + g λ µ ǫ νρστ = 0 (C. 24) contracted with the tensor P 1ν ǫ ρ (k 1 )ǫ σ (k 2 )k 2τ k 1λ followed by separate contractions with each independent momenta P 1µ , k 1µ and k 2µ we obtain three relations This means that only two tensor structures involving tensor ǫ αβγδ in (C.23) are independent and we choose as the independent tensors ǫ νk 1 k 2 P 1 and ǫ µk 1 k 2 P 1 . Thus the eq. (C.23) expressed in terms of these independent tensors has the form where again the scalar functions M

(i)
A depend on different Mandelstam invariants. Imposing the gauge invariant conditions k 1µ ǫ ν (k 2 )M µν A = 0 = ǫ µ (k 1 )k 2ν M µν A we obtain two relations Thus we can represent ǫ µ (k 1 )ǫ ν (k 2 )M µν A in two equivalent forms and A . (C.32) These expressions leads to the two following equivalent forms for the square M µν and M µν A (M Aµν ) * = 4(k 1 · P 1 ) 2 (k 2 · P 1 ) 2 |M (4) One can easily verify that similarly like in the case of decomposition of the vector tensor M µν V discussed earlier the expression (C.33) is directly obtained when one choses ǫ(k)·P 1 = 0 gauge in the eq. (C.28), whereas the expressions (C.34) correspond to the gauge choice ǫ(k) · k 1 = 0 in the eq. (C.28).

C.3 Relation with the present study
For our present study, the general results of the two previous subsections can be applied after the identification The P 1 -gauge of subsection C.1.3 is the axial gauge p · ǫ = 0 which was used in ref. [13], since in the collinear limit ∆ ∼ p.
D. Integration over z and x D.1 Building block integrals for the numerical integration over x Using the same notation as in ref. [13], we list for completeness the building block integrals which are involved in the numerical evaluation of the scattering amplitudes. For a generic GPD f, we define f (x, ξ) dx , (D.10) This set of 12 integrals is not minimal, and can be further reduced in terms of the 6 elementary integral I b , I c , I d , I e , I h , I i as follows. First, The other integrals simplifies when specifying the symmetry of the GPD f.
while for an antisymmetric GPD, one has (D.20) Each of the 6 elementary integral I b , I c , I d , I e , I h , I i is finite and is evaluated numerically, using our models for the various involved GPDs. After computing this set of integrals, the evaluation of the gauge invariant blocks of diagrams is straightforward using the decomposition given in two next subsections. Below, we will not indicate the function f , since it is obvious from the context.

D.2 Integration of gauge invariant sets of diagrams
We now present the result for the contributions of the various gauge invariant blocks of diagrams of figure 2 in terms of the 5 elementary integrals I b , I c , I e , I h , I i after integration over z and integration over x when multiplied by GPDs, which we denote generically as f q . One should note that the integral I d which appears in several diagrams, does not appear when considering gauge invariant sets of diagrams.

D.2.1 P P part
We decompose the trace involved in a diagram diag as where a prefactor C π as well as any charge coefficient has been factorized out. We denote the dimensionless coefficients This is in accordance to the conventions (5.10, 5.11) and (5.6, 5.7). For the block (AB) 123 made of diagrams For the block (AB) 45 made of diagrams These sums can be simplified when acting on GPDs with definite symmetries. For a symmetric GPD, using eqs. (D.14) and (D.15), we get

D.2.2 SP part
We decompose the trace involved in a diagram diag, as where a prefactor C π as well as any charge coefficient has been factorized out. We denote the dimensionless coefficients These definition are in accordance to the conventions (5.8, 5.9) and (5.4, 5.5).
For the block (AB) 123 made of diagrams For the block (AB) 45 made of diagrams These sums can be simplified when acting on GPDs with definite symmetries. For a symmetric GPD, using eqs. (D.14) and (D.15), we get and and Note that only the 6 coefficients       Let us discuss these various cuts with some details. First, since we rely on factorization at large angle, we enforce the two constraints −u ′ > (−u ′ ) min , and −t ′ > (−t ′ ) min , and take (−u ′ ) min = (−t ′ ) min = 1 GeV 2 . The first constraint is the red line in figures 11 and 12, while the second, using the relation M 2 γπ + t ′ + u ′ = t + m 2 π , is given by

PSfrag replacements
and shown as a blue line. The variable (−t) varies from (−t) min , determined by kinematics, up to a maximal value (−t) max which we fix to be (−t) max = 0.5 GeV 2 , these two boundaries being shown in green in figure 12.
The value of (−t) min is given by eq. (A.8). In the domain of M 2 γπ for which the phasespace is a triangle, as illustrated in figure 11, the minimal value of −t is actually above (−t) min . For a given value of M 2 γπ , this minimal value of −t is given, using eq. (E.1), by This constraint on −t leads to a minimal value of M 2 γπ , denoted as M 2 γπ crit , when (−t) inf = (−t) max , which thus reads With our chosen values of (−u ′ ) min , (−t ′ ) min and (−t ′ ) max we have M 2 γπ crit ≃ 1.52 GeV 2 , below which the phase-space is empty.
For the purpose of integration, we define, for −(u ′ ) min −u ′ , The maximal value of −u ′ , attained when −t = (−t) max , is given by figure 11.
With our choice of parameters, we get M 2 γπ trans ≃ 2.01 GeV 2 in the case of S γN = 20 GeV 2 .  Above this value, the phase-space is a trapezoid, illustrated in figure 12. This trapezoid reduces to an empty domain when (−t) min = (−t) max . This occurs for With our choice of parameters, we get M 2 γπ Max ≃ 9.47 GeV 2 in the case of S γN = 20 GeV 2 , a value which decreases with decreasing values of S γN .
The minimal value of S γN is obtained from the constraint M 2 γπ crit = M 2 γπ Max and equals S γN crit ≃ 4.75 GeV 2 .

E.2 Method for the phase space integration
The phase space integration goes along the same line as in ref. [13]. Using the phase-space described in the previous subsection, the integrated cross section reads dσ dM 2 γπ = θ(M 2 γπ crit < M 2 γπ < M 2 γπ trans ) (E.10) Using our explicit dipole ansatz for F (t), see eq. (6.2), we obtain , which is our building formula for the numerical evaluation of integrated cross sections.

F. Angular cut over the outgoing photon
In order to take into account limitations of detection of the produced photon, we compute the photon scattering angle θ in the rest frame of the nucleon target, with respect to the −z axis. We refer to appendix C of ref. [13] for details and only provide here the main formulas. At fixed value of M 2 γπ , we formally write tan θ = f (−u ′ ) . for tan θ < 0, θ = π + arctan(tan θ) .

(F.3)
For simplicity, we restrict our analysis to the ∆ t = 0 case.      The obtained angular distribution is shown in figure 13 for γπ + pair production on a proton target, and in figure 14 for γπ − pair production on a neutron target. These two angular distributions have a rather similar shape, with a clear dominance of moderate values of θ.  In practice, at JLab, in Hall B, the outgoing photon could be detected with an angle between 5 • and 35 • from the incoming beam.
The effect of an upper angular cut can be seen in figure 15 for γπ + pair production on a proton target, and in figure 16 for γπ − pair production on a neutron target. This effect is almost identical when comparing the γπ + and γπ − cases. As seen from figures 13 and 14, it mainly affects the low S γN domain. In particular, the effect of the JLab 35 • upper cut remains negligible as shown in figures 15 and 16, both for the γπ + and γπ − cases.