Perturbative Calculations of Gravitational Form Factors at Large Momentum Transfer

We perform a perturbative QCD analysis of the gravitational form factors (GFFs) of nucleon at large momentum transfer. We derive the explicit factorization formula of the GFFs in terms of twist-3 and twist-4 light-cone distribution amplitudes of nucleon. Power behaviors for these GFFs are obtained from the leading order calculations. Numeric results of the quark and gluon contributions to various GFFs are presented with model assumptions for the distribution amplitudes in the literature. We also present the perturbative calculations of the scalar form factor $\langle P'|F^2| P\rangle$ for pion and proton at large momentum transfer.


I. INTRODUCTION
The physics of the gravitational form factors (GFFs) has attracted a great interest of hadron physics community in recent years. The experimental investigations of the GFFs are one of the main goals for the ongoing JLab 12 GeV program [1] and the future electron-ion colliders [2][3][4]. The GFFs are defined as the form factors of the energy-momentum tensor (EMT), which describe the elastic interaction between the graviton and the particles [5][6][7].
The EMT of a field theory in the Minkowski space is derived from the space-time translation invariance of the action by the Noether theorem [8]. The EMT obtained in this way is called the canonical EMT and normally not symmetric with respect to its Lorentz indices. However, the physical concept of the GFFs actually stems from another kind of the EMT which plays a unique role in the Einstein theory of gravity [9]: where g µν (x) is the classical gravitational field in the general relativity and S M is the matter action of the system. This form of EMT is automatically symmetrized and differs from the canonical one by a total derivative, hence it is called symmetric EMT or Belinfante-improved EMT [10,11].
To explain the implication of GFFs, we can consider a hadron system in the classical gravitational background. When the curvature of the gravitational field is weak enough, we can expand the metric field around the Minkowski metric η µν = (1, −1, −1, −1), Then what we obtain is a flat-space action in which the perturbation field h µν (x) is equivalent to a spin-2 field and the corresponding particle can be considered as the graviton. In the action, it is interesting to find that the graviton couples with the EMT operator in the following form: That means that when a hadron N scatters with the graviton elastically, the only way the quarks and gluons inside the hadron can interact with them is through EMT operator. Therefore, this scattering is describe by the following transition matrix N (P )|T µν QCD |N (P ) .
Although it is not realistic to probe the GFFs directly through the graviton-hadron scattering in experiments, they can be alternatively extracted from the generalized parton distributions (GPDs) in the hard exclusive processes [7,[38][39][40][41] , e.g., from the deeply virtual Compton scattering [7,38,[42][43][44] and the deeply virtual meson production [45][46][47]. Meanwhile, they can also be explored in the time-like processes such as the hadron productions in the two-photon collisions, see, for example, a recent analysis [48] on the experimental data from the Belle Collaboration at KEK [49].
From the theory side, the GFFs can be computed nonperturbatively in the framework of the lattice QCD, see, e.g., the pioneer works in Refs. [28,[65][66][67]. On the other hand, if the momentum transfer (−t) is sufficiently large, one can also carry out the perturbative analysis thanks to the asymptotic freedom of the QCD. This is the main subject of this paper. We will perform a systematic analysis on the complete set of the nucleon GFFs at large momentum transfer, including both quark and gluon sectors. The methodology is based on the widely-used Efremov-Radyushkin-Brodsky-Lepage (ERBL) formalism in hard exclusive processes [68][69][70][71][72][73][74] and the progresses in the classification of the nucleon light-front wave function involving the orbital angular momentum [74][75][76] as well as the power counting rule [77][78][79]. We will mainly focus on the nucleon GFFs and derive their factorization formulas in terms of twist-3 or twist-4 nucleon light-cone distribution amplitudes [80,81]. Parts of our study on the gluonic GFFs as well as the pion GFFs have been reported in the letter recently [82].
The rest of the paper is organized as follows. In Sec. II, an introduction of the quark and gluon EMT and the GFFs will be given. In Sec. III, we perform the perturbative analysis of the nucleon GFFs at large momentum transfer. In Sec. IV, we will present the numeric evaluation of the individual GFFs based on model assumptions for the twist-3 and twist-4 distribution amplitudes. In Sec. V, the scalar form factors at large (−t) will be discussed for both pion and proton. Sec. VI is the conclusion of the paper.

II. GRAVITATIONAL FORM FACTORS
In this section, we will first present the explicit expressions of the symmetric QCD EMT and discuss some important features of it.

A. Energy-Momentum Tensor
The symmetric QCD energy momentum tensor is defined as The quark and gluon sectors are given by where ψ f is a quark field of flavor f , and F a µν = ∂ µ A a ν − ∂ ν A a µ − g s f abc A b µ A c ν is the strength tensor of the gluon field A a,µ , and the covariant derivatives are defined as The sum of the quark flavors and colors are implied.
There are several features that we would like to emphasize. First, the total EMT in Eq. (5) is conserved due to the space-time translation invariance, i.e. ∂ ν T µν QCD = 0. However, the gluon or quark part of the EMT is not conserved individually. Instead, the derivative of the EMT operators are related to the twist-four operators [83,84], The above equations are derived by applying the equations of motion for the quark and gluon fields: With the EOM of the gluon strength tensor, one can also find the right sides of Eqs. (8) and (9) cancel out each other. This confirms the conservation of the total EMT operator explicitly. Therefore, the total EMT operator is UV finite and scale-independent [85][86][87][88] while the quark or gluon sector of the EMT is divergent and needs regularization and renormalization [17,18]. The renormalization of the EMT is closely related to the trace of this operator. Classically, the quark or gluon EMT is traceless if the quark masses are neglected: which is also the consequence of the conformal symmetry in the massless QCD. This feature will be broken at the quantum loop level, called trace anomaly [13,[85][86][87][88][89]: where β(g s ) is the QCD beta function and the terms (1 + γ m )m fψf ψ f should be added in the massive case.
Second, in the covariant quantization of QCD, the EMT operators in Eqs. (5,6) are not completed rigorously. For example, in the Faddeev-Popov quantization procedure for the gauge theory one has to introduce the ghost and gauge-fixing terms in the Lagrangian, thus there are corresponding terms in the QCD EMT [85]. However, these terms are ensured to be absent in the physical matrix elements by the BRST symmetry [90] and hence have no observable effects. Therefore, we have excluded them in our study of the GFFs.

B. Nucleon Gravitational Form Factors
Nucleon is a spin-1/2 particle, and the quark or gluon gravitational form factors can be parameterized as [7]: where a = q, g and the brace for the Lorentz indices µ, ν denote the symmetrization known as A (µ B ν) ≡ (A µ B ν + A ν B µ )/2. P and P are the momenta of the initial and final particle respectively, satisfying the onshell condition P 2 = P 2 = M 2 , where M is the nucleon mass.P is the average momentum, ∆ is the momentum transfer and t is the momentum transfer squared: The nucleon state is covariant normalized, P |P = 2P 0 (2π) 3 δ 3 ( P − P ). u s (P ) is the Dirac spinor of nucleon normalized asū s (P )u s (P ) = 2M δ ss . One can utilize the well-known Gorden identity for the spinor, to obtain another parametrization for the GFFs: where J a (t) = [A a (t) + B a (t)]/2. Here we follow the notations in Refs. [7,38] for the C form factors. They have also been referred as D or d 1 form factors in Refs. [26-28, 33, 34] with different normalizations: From the above expressions, one can observe the quark and gluon GFFs are the Lorentz invariants of the momentum transfer squared t. As shown in Eqs. (8,9), the quark and gluon parts of EMT are not conserved and require renormalization, and the corresponding GFFs should also depend on the renormalization scale [17,18]. For each form factor, we can also define the total GFF, e.g. C(t) = a=q,g C a (t). They are, on the other hand, renormalization scale independent due to the conservation of the total EMT. This feature has further implication on the GFFs: where the Hesienberg equation of the EMT operator is applied. The last equality yields where the terms related to the A, B, C-form factors vanish automatically from its own tensor structure or the on-shell condition of the Dirac spinor, / P u(P ) = M u(P ). Therefore, an important constraint on the C-GFFs is obtained: This constraint works regardless of the value of the momentum transfer squared t and will serve as an important consistent check on our perturbative calculations of the quark and gluon GFFs at large momentum transfer. The details will be presented in Sec. III F. The GFFs C q,g (t) play an important role in the nucleon mass reconstruction since it is closely related to the trace anomaly F 2 and its value at t = 0 determines the contributions of the so-called quantum anomalous energy [12-15, 17, 55]. The C q,g (t) form factors in Eqs. (12,15) are proposed to describe the mechanical properties such as pressure and shear forces distribution inside the nucleon [26,27,31,32]. This interpretation uniquely requires the whole knowledge of t-dependence from C-form factor.
In addition, the A q,g form factor at zero momentum transfer can be interpreted as the longitudinal momentum fraction that the parton carries inside the nucleon in the infinite momentum frame with a=q,g A a (0) = 1. The J-form factor J a (t) = [A a (t) + B a (t)]/2 at t = 0 describe the partitions of the quark and gluon angular momentum from the nucleon spin satisfying the constraint a=q,g J a (0) = 1 2 , which is known as the Ji's sum rule in the literature [7].
In the following sections, we will carry out a systematic investigation on the nucleon GFFs for the quarks or gluons at large (−t).

III. PERTURBATIVE ANALYSIS AT LARGE MOMENTUM TRANSFER
In general, the nucleon GFFs are non-perturbative due to the color confinement. Nonetheless, in the large momentum transfer limit, there are indeed perturbative calculable effects in the GFFs. The physics behind is that the highly virtual graviton has short enough wavelength to resolve the structure of the nucleon by the short spacetime interactions. The partons participated in this interaction are weakly coupled but endure with hard gluon exchanges to pass the large momentum transfer from the initial nucleon to the final nucleon. These effects can be calculated using perturbation theory due to the asymptotic freedom of QCD. On the other hand, the interac-tions among those active partons and the other spectator partons inside the nucleon are still of long range and hence strong. In this section, we will demonstrate that this short-range and long-range physics can be factorized at lowest order perturbation theory, leading to the factorization theorem of the GFFs in terms of twist-3 or twist-4 nucleon light-cone distribution amplitudes. Our derivations follow the ERBL formalism [68][69][70][71][72][73][74] with development on the nucleon light-front wave function involving zero or one unit of orbital angular momentum [74][75][76].
We will first introduce the helicity-amplitude with appropriate tensor projection to isolate each GFF, and choose a reference frame. Later, we will review the threequark Fock state of the nucleon in the literature. Based on these materials, we will carry out the perturbative analysis of the GFFs through the helicity-conserved and helicity-flip amplitudes, respectively. Finally, we summarize the factorization results and perform consistent checks.

A. Helicity Amplitude of EMT
Since the nucleon is a spin-1/2 particle, the parametrization of the EMT amplitude involves the Dirac structures as well as the tensor structures. The corresponding form factors can be extracted from the nucleon helicity amplitude with appropriate tensor projections. In the high energy, we will find that the helicity-conserved amplitudes yield the A-form factors, and the helicity-flip amplitudes lead to other GFFs.
As presented in Eq. (15), there are two Dirac structures involved in the EMT amplitude:ū s γ µ u s andū s u s , corresponding to the chiral-odd and chiral-even Dirac bilinear spinors. In the large momentum transfer limit, the nucleon mass can be neglected, and hence the chirality becomes a good quantum number to identify the nucleon state, i.e., where u ↑/↓ (P ) is the chirality eigenstates, named as right-or left-handed spinors. Besides, for a massless spin-1/2 particle, its chirality is known to be equivalent to its helicity, where helicity is defined by the projection of the spin along the three-momentum direction. We will use helicity as the terminology in the following paragraphs. Now we choose the appropriate helicity eigenstates for the initial and final nucleons along with tensor projections to extract the GFFs. Helicity-conserved Amplitude A-form factor can be extracted form the helicity-conserved amplitude: where B-form factor are power-suppressed as shown later and hence neglected in the above equation. Since only A-form factor is relevant at the leading power, we do not need to assign the tensor projection. Here we used the amplitude with positive helicity. The amplitude for the negative helicity can also be used and leads to the same result.
Helicity-flip Amplitude On the other hand, helicityflip amplitude can be used to extract the GFFs B a , C a and C a : To further isolate the GFFs, we can define the following projection tensors, Notice that g µν Γ µν 1,2 = 0 in the massless limit. Applying the above, we can derive From the traceless feature of the EMT in our calculations at this order, we can also derive C form factor: As mentioned above, due to the conservation of the EMT, the total C-form factor should vanish. That renders Furthermore, the total contributions to the B and C form factors follow As a consequence, we find that there is a nontrivial relation between the total B-and C-form factors at large momentum transfer at the leading order perturbation theory: This relation along with the constraint C(t) = 0 will serve as the consistent checks of our calculations.

B. Breit Frame
In order to facilitate the computation, we need to fix the reference frame explicitly. Since the GFFs are the Lorentz invariant, one can choose any frame in principle. Particularly, we use the so-called Breit frame in the calculation where the initial nucleon is along the z-axes with high energy and the final nucleon is along the opposite direction: Here the large momentum transfer limit −t M 2 is applied and the light-cone coordinates are used where In the Breit frame, the momenta of the initial state and final state nucleons as well as the momentum transfer contain only the longitudinal components. Hence, all the large scales in the partonic subprocesses must come from the longitudinal momentum and make the power counting of the perturbative part more transparent.

C. Three-valence Quark Fock State of the Nucleon
In the Breit frame, both the initial and final nucleons are traveling along the longitudinal direction and carrying large momentum. In this frame, one can apply the light-front Fock state expansion [91]. In general, there are infinite number of Fock states with light-front wave functions in this expansion. However, for an exclusive process at the large momentum transfer, one can conclude from the general power counting rule that the leading contributions come from the terms with minimal numbers of partons and the fewest orbital angular momentum (OAM) content [77][78][79]. In the following, we will review the three-valence quark Fock state of the nucleon [75,76].
For the nucleon, the leading Fock state is made of three valence quarks, i.e. u, u, d. Each quark can carries the helicity λ i = ±1/2. The total helicity of the valence quarks has the value as 3/2, 1/2, −1/2, −3/2 and the corresponding Fock states can be denoted where the upper subscript represents the total quark helicity λ = i λ i and the superscript stands for the projection of OAM along the z direction determined from the angular momentum sum rule l z = Λ − λ. Similarly, the left-handed proton can be expressed as Therefore, the proton helicity amplitude of EMT can be expressed with the above Fock state components. Then the nucleon EMT amplitudes are factorized into the partonic EMT amplitudes multiplied by the lightfront wave function from the initial and final state nucleons. In the high energy scattering, the light quark mass can be neglected, and the parton helicity is conserved. That means the nucleon helicity-flip amplitude can only be obtained by involving a non-zero OAM from either initial or final state Fock components. In fact, this OAM along the longitudinal direction is achieved by the parton's intrinsic transverse motion and hence introduces some transverse momentum factors in the phase space integral which normally scale as Λ QCD . In the end, this factor will be compensated by the large momentum scale (−t in our case) and lead to power suppression. Therefore, more OAM involved, more power suppression will be obtained [79].
For the nucleon helicity-conserved amplitude in Eq.(20), we have where the power suppressed contributions from ∆l z = l z −l z = 0 are neglected. However, for the nucleon helicity flip amplitude in Eq. (21), one unit of OAM is needed to achieve the spin flip. Hence, the leading contributions come from the components with ∆l z = 1, where the helicity flip for the first term is from the initial state OAM and the second one from the final state OAM. From the parity and time reversal invariance symmetry, one can show that they in fact have the same contributions to the GFFs at large −t. Here we only focus on the first term.
In the followings, the expressions of the relevant Fock state components involved in our evaluation for the nucleon helicity amplitude are given [75,76]: where the three-valence-quark state is defined by In the above expressions,û † a,λ andd † a,λ are the creation operators of the u, d quarks with the helicity λ and color indices a in the fundamental representation of SU(3) group, normalized as {û a,λ (k),û † b,λ (k )} = δ λλ δ ab 2k + δ(k + − k + )δ (2) (k − k). ψ 1 , ψ 3 , ψ 4 are the light-front wave functions of the proton depending on κ i ≡ (x i , k i ): x i is the longitudinal momentum fraction of the parton with i x i = 1; k i is the intrinsic transverse momentum of the parton with i k i = 0. The integral measures for the phase space are defined as The momentum factor k L ≡ k x − ik y is the explicit representation of quark OAM in the proton.

D. A-Form Factor and Helicity-conserved
Amplitude at Large −t As introduced in the Sec. III A, the A-GFFs are related to the nucleon helicity-conserved amplitude. The dominated contributions come from three-quark light-cone wave function with zero orbital angular momentum: At large momentum transfer, the transverse momenta of the partons in the above amplitude can be neglected. Therefore, where Φ 3 is the twist-three light-cone amplitude of proton [80]: The label uud and udd denote the relevant Fock state configurations in Eq. (39): Since the helicity and flavor of a massless quark are conserved in the fermion line, there is no interaction between these two configurations. Both of them share the same kind of diagrams without any quark line crossing, as presented in Fig. 1. The second configuration in Eq. (44) have additional diagrams with the final u quarks exchanged from Fig. 1. Therefore, the hard partonic part for the EMT amplitude can be expressed as where and H a is obtained from H a by interchanging y 1 and y 3 . As shown in Fig. 1, for the quark contribution, the numbered places represent the insertion of the quark EMT operator. Similarly, we can have gluon contribution, where the insertion only appears on the two gluon propagators. The total contribution to the GFF comes from the sum of all these diagrams with all possible insertions. The diagrams in Fig. 1 represent the two different configurations: (1) two gluons attach to the spin-down quark; (2) two gluons attach to the spin-up quark. All other diagrams can be obtained by an exchange of the quark lines (1 ↔ 3).
Since the partons are approximately on-shell, we use the covariant perturbation theory in Feynman gauge to compute the partonic amplitudes. The technical part of the computation is to evaluate the Dirac structure which involve three fermion lines. In general, they have the following forms:ū Since we only concentrate on the Dirac structure, all irrelevant factors are included in the Γ 1,2,3 , which denote the product of γ-matrix with odd number from three fermion lines presented in the Fig. 1. All the diagrams share the same color factor, where T a is the generator in the fundamental representation. Note that the helicity amplitudes in Eq. (47) have a pair of quark lines with zero total helicity. One can apply the following strategy to combine them into a Dirac trace. First, we move out the momentum faction inside the spinor by u λ (x i P ) = √ x i u λ (P ), and similarly for u λ (y i P ). The next step is to utilize the identitȳ where Γ R is obtained from Γ by reversing the order of the γ-matrix. Then we can use the spin density matrix for massless spinor which will turn two fermion lines with opposite helicity of Eq. (47) into a Dirac trace. After that, Eq. (47) becomes Evaluating the trace, contracting the Lorentz indices, and applying the on-shell condition for the spinor / P u λ (P ) = 0,ū λ (P ) / P = 0, we arrive at the following four Dirac bilinear: where In particular, from the explicit calculation, one can find that E 3 and E 4 have the form as x i − y i , and hence the second line with ∆ factor should vanish due to the symmetry between the initial and final state. Furthermore, one can eliminate the Levi-Civita tensor by the identity: which can be verified by using the explicit expression of the Dirac spinors given in Eq. (A4). Therefore, one finally obtain the Dirac structure in the proton helicityconserved amplitude as expected: The hard coefficients for A a (t) are straightforward to find out. After summarizing the results, the A a (t) GFF for the nucleon has the following factorization formula at large momentum transfer: is the twist-three light-cone amplitude of the proton [80]. The hard function can be written as where A a is obtained from A a by interchanging y 1 and y 3 . C B = 2/3 is the color factor. For the gluon part, A g can be written as For the quark part, the hard coefficient can be expressed as In summary, the quark and gluon A-form factors have the same power behavior at large momentum transfer.
E. Helicity-flip Amplitude at Large −t Now we turn to the computation of B, C and C form factors at large −t. As shown in Eq. (25), they can be extracted respectively from the proton helicity-flip amplitude of EMT with appropriate tensor projections. The dominant contributions come from the amplitude components with one unit of OAM from either initial or final state nucleons. Since they yield the same contribution, we only focus on one of them: where the factors k L = k x − ik y are manifestations of the quark OAM with l z = −1 from the initial state proton. The partonic amplitudes have the same structure as that in Eq. (45) except that the transverse momenta of the initial state partons can not be ignored. Otherwise, the amplitudes will vanish due to the OAM factor k L in the phase-space integral. Instead, we need to evaluate the following partonic amplitude: where we have used the label Γ µν to represent the tensor projection for the B, C, or C form factors in Eq. (25). The partonic state |uud has been defined in Eq. (44). Similar to the A-form factor, there are three fermion lines responsible for the three-valence-quark configuration for each diagram in Fig. 1. These perturbative diagrams share the same color factor as that in Eq. (48). They also contain the similar Dirac structure as that in Eq. (47) so that we apply the similar strategy for the Dirac algebra, following Eqs. (49)(50)(51) and Eq. (53). In order to obtain the leading contribution at large −t, we should expand the transverse momenta in terms of k i / √ −t. The linear terms of k i will contribute at this order. In particular, one can apply the following formula for the spinor: which leads to a linear expansion term. Applying the transverse momentum expansion of the partonic scattering, we obtain the following expression for the linear expansion term, where k 2 = −k 1 − k 3 has been used to simplify the final results. To obtain the helicity-flip like structure, one can use the identity γ i u ↑ (p) = (δ ix + iδ iy )u ↓ (p) for i = x, y, and hence the bilinear spinor can be further reduced as u ↑ (P )/ k i u ↑ (P ) = k iRū↑ (P )u ↓ (P ). Therefore, the GFFs can be expressed as where k R/L = k x ± ik y . The hard parts have the form as where the x i , y idependences from the phase integral are also included. On the other hand, from the property of the transversemomentum integral, one has [ The transverse-momentum dependence from the initial state can be absorbed into the twist-4 light-cone distribution amplitudes of nucleon by the following relations [74]: The normalizations of these twist-4 amplitudes are the same as those in [81]. Similarly, the transverse momenta of the final wave function can also be integrated out as in Eq. (42) and thus the twist-3 distribution amplitude of the final nucleon is obtained. Therefore, the GFFs at large momentum transfer have the following general factorization structure: where Φ 3 is the twist-3 distribution amplitude and Φ 4 , Ψ 4 are the twist-4 distribution amplitudes defined above. They embrace the non-perturbative effects of the GFFs at large −t. H Φa , H Ψa are the hard coefficients and can be calculated order by order in principle.
Although we derived the above formula at the lowest order perturbation theory, we expect that it still takes the same form when high-order α s -corrections are considered. We notice that the bare quark or gluon GFFs suffer from UV-divergence and need renormalization at higher orders [17,18]. The operator mixing is needed and the trace anomaly F 2 plays an important role in the renormalization of EMT operator [13,17,18,[85][86][87][88][89], which is not present in our leading order calculation. As shown in Sec. V, the P |F 2 |P for nucleon at large −t is also related to the helicity-flip amplitude. Hence, the factorization structure of B, C, C form factors in Eq. (65) should not be altered by the operator mixing.

F. C Form Factors and Check on the Cancellation between Quarks and Gluons
In the last section, we have explained the general procedure on the calculation of the B, C and C form factors from the helicity-flip amplitude with different tensor projections. In particular, C form factor is actually determined by the traceless feature g µν T µν q,g = 0 at this order. It is an important cross check that the total contribution of C from the quarks and gluons vanishes because of the current conservation. In order to show that from the explicit calculation, we define the following partonic amplitude for C-GFFs: where Γ µν 2 = ∆ µ ∆ ν − ∆ 2 g µν /4 has been defined in Eq. (25).
Following the strategy presented in the last section, we can go through the detailed but straightforward computation on the perturbative diagrams in Fig. 1 and check explicitly the cancellation of C form factor contributions from the quarks and gluons. In the followings, we will list the contributions of C q,g and an overall factor is implied for brevity: First, we summarize the contribution to the gluon GFF, The total contribution will be obtained by adding the above together and apply (1 ↔ 3).
To compute the quark contributions from the diagrams of Fig. 1, we notice that the total contribution from g µν term in Γ µν 2 cancel out among the diagrams. Therefore, we only need to take into account the ∆ µ ∆ ν term in Γ µν 2 in the projection of these diagrams. For the first diagram, we have It is interesting to find out that the total contribution from the quark and gluons takes the following expression, which is anti-symmetric under the exchange of (1 ↔ 3). Therefore, the cancellation of C form factor between the quark and gluon contributions happens when we adding the first diagram of Fig. 1 and its mirror. For the second diagram in Fig. 1, we have We can also sum up quark and gluon contributions in the second diagram, which does not vanish, even considering (1 ↔ 3) symmetry. However, we will show this is canceled by the contributions from the third diagram. The quark contributions from the third diagram of Fig. 1 are very simple, and all others vanish. Adding them together, we have, This exactly cancels that from C q,g . To summarize, we have shown that the total C form factor factor vanishes when adding all contributions from the quarks and gluons. This provides an important cross check of our derivations. Moreover, from above results, we can obtain the factorization formula of C a form factor for quark or gluon respectively. It takes the same form as shown in the Eq. (65): where the hard functions have the following structure: and C is obtained from C by interchanging y 1 and y 3 . At the leading order perturbation theory, we obtain The functions K i are defined as

G. B and C Form Factors at Large −t
Similar to the C-form factor, the B and C-form factors for quarks and gluons come from the helicity-flip amplitude and can be calculated by the same analysis explained in Sec. III E at large momentum transfer. They all follow the same factorization structure as in the Eq. (65) or Eq. (75). In this subsection, we will summarize these results.
a. Gluon form factor With the above analysis, we carry out a detailed derivation for all the diagrams of Fig. 1 and C g (t) can be factorized into, The hard functions can be written as where C is obtained from C by interchanging y 1 and y 3 . From the detailed calculations of the diagrams in Fig. 1, we obtain where and functions K i defined in Eq. (78).
Likewise, the gluonic B-form factor has the similar factorization formula: where the hard function H for B g is and B is obtained from B by interchanging y 1 and y 3 . Detailed calculation gives at the leading order.
b. Quark form factor We turn to the case of the quark form factors. Similar analysis and derivation show that C q (t) can be factorized as where Ψ 4 and Φ 4 are the twist-four distribution amplitude of the proton [81]. The hard fucntions have the structure: where C is obtained from C by interchanging y 1 and y 3 . From the detailed computation of diagrams in Fig. 1, we obtain Similarly, the B-form factor for quarks follow the same form of factorization formula at large −t as the gluonic one and the corresponding hard fucntions have the same structure. Explicit calculation yields The functions H qC and H qB can be unifiedly described by the function H q with different signs chosen. The upper sign is for the C q , while the lower sign is for the B q , i.e. where The functions G i are defined as There are several points we would like to comment on the above results. First, one can carry out the consistent check by calculating the C a from factor We have checked that the obtained results are the same as those in Eq. (75,76,77) and hence the total C GFF vanishes as expected, see also, the relation between the total B and C form factors of Eq. (30). Second, the above results confirm the power counting for the GFFs: C q,g , B q,g ∼ 1/(−t) 3 , C q,g ∼ 1/(−t) 2 . Compared with the A-form factor, the B, C-form factor are 1/(−t)-suppressed due to the quark OAM for the helicity flip amplitude. C q,g also come from the helicity-flip amplitude, but with one power of (−t) enhanced compared to B q,g , C q,g . This is because the parameterizations of the GFFs.
Third, it is interesting to see that the perturbative coefficients for B q,g , C q,g form factor can be described by one function H q,g for quark and gluon, respectively: B q,g and C q,g can be related to each other only by reversing the sign of some terms. In addition, the hard coefficients of B q,g , C q,g , C q,g GFFs suffer from the end-point singularity when x i , y i → 0. This effect comes from the soft partons carrying small fraction of the longitudinal momentum. Since the longitudinal and transverse momenta of the partons are comparable, we are not able to perform the intrinsic transverse momentum expansion in the end point region. One should develop a rigorous framework to factorize and resum these soft parton contributions in the GFFs, which however is beyond the scope of the current paper. Nonetheless, we would like to point out some observations on the end point behaviors from our results. From Eqs. (81,82,85), we find that the end point singularities of the gluonic B and C GFFs can be related by and hence (95) from Eq. (93). This relation dose not hold for the quark contributions and the end point singularity does not cancel in the total quark and gluon contributions. However, from C q = −C g , we also have One more point, the GFFs can be obtained from the GPDs H a , E a by taking the first moment. For the gluonic GFFs, we have carried out a consistent check by using the gluon GPDs at large −t presented in [59]. For the A q GFF, we use the quark GPDs H q in [92] to check the result, whereas for the other quark GFFs, we compute the quark GPD E q by using the methods introduced in the previous sections. We list them in the Appendix for reference. The results are all consistent.

IV. NUMERICAL RESULTS
In this section, we will present numeric results for the gluon and quark GFFs of the nucleon at large (−t). To achieve that, the non-perturbative inputs of the twistthree and twist-four nucleon light-cone distribution amplitudes (DAs) are needed [80,81]. There have been progresses in the determination of the nucleon DAs with different methods in the literature, for example, QCD sum rule [81,[93][94][95][96][97][98], lattice simulation [99][100][101][102][103][104], and the phenomenological models [95-98, 105, 106]. Several works  fN , λ1, ϕ10, ϕ11, η10 and η11. The values are presented at the renormalization scale µ 2 = 1GeV 2 by using one-loop evolution. The errors are not shown here. The values of asymptotic (Asy.) model and QCD sum-rule (SR) estimate are taken from [98]. Especially, it is well-known that the parameters here are related to the local quark-gluon operators which the QCD sum-rule methods can be applied in [81, 93-98, 111, 112]. The Braun-Lenz-Wittman (BLW) model is also shown here [98]. This model is based on the light-cone sum-rule analysis on the nucleon electromagnetic form factors which results in a relation between the electromagnetic form factors and the nucleon twist-3 and twist-4 DAs. With the choice of specific values, a good description for the experimental data was obtained. Furthermore, the parameters have been estimated in the lattice QCD simulations [99][100][101][102][103][104] and the updated values evaluated by RQCD collaboration are adopted [104].
In our calculations, we adopt the parametrization of the nucleon DAs proposed in [80,113]. For the twist-3 DA, we use . For the twist-4 DAs, we use where W Φ,Ψ are known as the the Wandzura-Wilczek contributions [113], which can be expressed in terms of the parameters presented in the twist-3 DA: while R Φ,Ψ contain new parameters: In the above equations, the scales of the parameter are hidden for brevity, e.g. f N ≡ f N (µ 2 0 ). The label L stands for the one-loop evolution factor from the scale µ 0 to µ: The running coupling at the leading order is where b = 11 − 2/3n f and n f is the number of the active quark flavors. Here we take n f = 4 and Λ QCD = 154 MeV. Other parameterizations are also proposed in [81,111,112]. To our accuracy, there are total of 6 non-perturbative parameters: f N , λ 1 , ϕ 10 , ϕ 11 , η 10 and η 11 . f N , λ 1 are the normalizations of the nucleon DAs [113]. ϕ 10 , ϕ 11 , η 10 and η 11 are the shape parameters which describe the shape of the amplitudes and also determine the deviation from the asymptotic forms. Those parameters can be estimated by using the non-perturbative methods [81,[93][94][95][96][97][98][99][100][101][102][103][104]. Phenomenological models have also been built from fitting to the experimental data [95-98, 105, 106]. For the detailed summary and comparison of different estimates, we refer the readers to [101,104,106,108,112] and the references therein. In Table I, we list those applied in this paper.
With the above parametrization, we can carry out the convolution integral in the factorization formula for the GFFs. The calculation is straightforward for the A q,g form factor. However, as we mentioned above, for the B q,g , C q,g , C q,g form factors from the helicity-flip amplitude, there are end-point singularities in the convolution integral at x i , y i → 0. To regulate these divergences, we introduce a cut-off Λ 2 c /(−t) for the integration as proposed in [74] for Pauli form factor, where Λ c is a soft scale. With the leading-order results derived here, one can check explicitly that the end-point behaviors of integrand can only follow the form: 1 xi , 1 yi or 1 xiyj . Contributions from the first two forms will give the single logarithm term as ln(−t/Λ 2 c ) and the last one would yield the double logarithm term ln 2 (−t/Λ 2 c ). With that improvement, we find that the B q,g , C q,g scale as ln 2 (−t/Λ 2 c )/(−t) 3 at large momentum transfer, same as their total GFFs. C q,g scales as ln 2 (−t/Λ 2 c )/(−t) 2 . Similar double logarithms were observed in the pQCD analysis of Pauli form factor F 2 and play an important role in describing the scaling behavior of Q 2 F 2 /F 1 at the large momentum transfer as compared to the experimental data [74,114].
In Fig. 2 we show the (−t)-dependence of the GFFs A q,g (t) and A(t) = A q (t) + A g (t) obtained from different estimates of twist-3 nucleon DAs. Several features can be found. First, it is easy to see that all A-form factors are positive, which show the same sign with the lattice simulations of GFFs [65,67] at low −t. Second, different models differ significantly at low (−t). However, at higher (−t), the QCD evolution tends to reduce the model dependence for the distribution amplitudes and all predictions are within the same order of magnitude. Third, among the models, the A-form factor from the asymptotic DA has the smallest magnitude over the whole range of −t, while the DAs obtained from the QCD sum-rule yield the largest magnitude. Moreover, the numerical evaluations show the following relations between the A-form factors: A(t) > A q (t) > A g (t). For example, in Fig. 3, we present the comparisons of A q,g (t) and A(t) with the DAs obtained from the QCD sum rule. Similar relations can be observed in the previous lattice calculations of GFFs in the range 0.3GeV 2 < −t < 2GeV 2 [65,67].
In Figs. 4-6, we show the numeric results for B, C and C form factors. These GFFs are extracted from the nucleon helicity-flip amplitudes. To regulate the end-point singularities, we introduce a soft scale Λ c . First, we can make the comparisons among the quark, gluon and total GFFs with a fixed Λ c . Some general features can be found for the models in Table I. For the C-form factors, C q (t) and C(t) are negative while C g (t) is positive, and they have the relations |C q (t)| > |C(t)| > |C g (t)|. On the other hand, the lattice calculation in [67] shows that C g (t) is negative at low −t. That means the C g (t) will change sign at higher (−t). For the B-form factors, they are all negative with the relations |B g (t)| < |B q (t)| < |B(t)|. For the C-form factors, the total Cform factor vanishes as we confirmed explicitly before and C q (t) = −C g (t) > 0. In particular, we display the above comparisons among the GFFs with the nucleon DAs from QCD sum rule [81,98] in Fig. 4 where Λ c is taken as 200MeV.
On the other hand, from our perturbative results at leading order, one can check explicitly that the combination B g (t) + 4C g (t) is indeed free from end-point sin-

V. SCALAR FORM FACTORS
The scalar form factors(FFs), defined as the hadron form factors of the F 2 operator, are also important in GFFs physics, since it has the close connection to the trace anomaly, as shown in Eq. (11). Following the discussion in the previous sections, similar perturbative analysis can be applied to the scalar FFs at large −t. We will first study the pion case, which will be simple and interesting, and then turn back to the proton case.

A. Pion
The scalar FF for pion, G π (t), is defined from the parameterization of the F 2 transition matrix: At large (−t), the factorization analysis can be applied. At the leading order of α s , there is only one diagram that contributes, see, for example, Ref. [82]. Therefore, where q 1 = xP − yP , q 2 =xP −ȳP are the momenta that flow into the F 2 operator. φ(x) is the pion twist-2 light-cone distribution amplitude, defined as with the normalization , where f π ≈ 92.3MeV is the pion decay constant.
From the above results, it is interesting to see that the pion scalar form factor has no power behavior with respect to −t, which is different from the pion A, C form factor with A π g,q ∼ C π g,q ∼ 1/(−t) [82]. It is also expected that higher order contributions will not change the power behavior. At leading order perturbation theory, the possible (−t)-dependence only comes from the running of the strong coupling α s with µ 2 taken as −t. There are also evolution effects hidden in the pion DAs. Hence, to make a sensible numerical estimate on the large −t behavior of the pion scalar form factor, more details on the pion DAs are needed. Similar to the proton DAs introduced in the previous sections, the pion DA can also be expressed as the series in the the eigenfunctions of the leading order evolution equation (e.g. [115]): where C α n (x) are the Gegenbauer polynomials. L n (µ 2 ) has been given in Eq. (103) and represents the logarithmic factor for the one-loop evolution from the default scale µ 2 0 to the varied scale µ 2 . γ n is the anomalous dimension given by where ψ(n) is the digamma function and b = 11 − 2/3n f . When µ 2 approaches to +∞, all the logarithmic factors are suppressed, and the first term determines the asymptotic form of the DA, i.e., φ asy (x) = √ 6f π x(1 − x). The coefficients a 2n , known as Gegenbauer moments, describe the deviation from the asymptotic DA.
Following Eq. (11), it is convenient to express the pion transition matrix element of the trace anomaly in terms of the pion scalar form factor in massless-quark limit. With the Gegenbauer expansion of the pion DAs in Eq. (106) and µ 2 = −t, a simple expression can be obtained where only the sum of the Gegenbauer moments is relevant. The above result shows that the pion trace anomaly matrix element should be negative in the large −t limit. However, it is known that at zero momentum transfer, the pion trace anomaly should be positive [12,22]. The above results reveal the fact that the trace anomaly matrix element should change sign from the large (−t) to small (−t). Preliminary lattice simulation also indicates this behavior [116]. Sign changes have also been observed in the gluon C form factor for both pion and proton. Later we will also show the proton scalar form factor experiences the sign change.
With a truncation on the Gegenbauer-polynomial expansions and the explicit setting on the values of the coefficients a 2n at the default scale µ 2 0 , several models on the pion DA can be obtained. Here we follow the models adopted in [117]. In Model I [118], the coefficients are set as where µ 0 = 1GeV and the uncertainties are neglected here. These values are fitted from the QCD light-conesum-rule results for the pion electromagnetic form factor to its measured data [118] . In model II [119][120][121], the coefficients are set as with µ 0 = 1GeV. The values are taken from an estimation from QCD sum rules with nonlocal condensate [119][120][121][122]. Besides the Gegenbauer-type models, the Ads/QCD inspired model(Model III) are also applied [117,123]: with α π (µ 0 ) ≈ 0.422.

Asy.
Model I

Model II
Model III Adopting the above models, we conduct the numerical evaluations on the pion gluon trace anomaly transition matrix. The results are displayed in Fig. 7 with respect to (−t) in the range from 1 GeV 2 to 5 GeV 2 . Here we take n f = 4 and Λ QCD = 154MeV.

B. Nucleon
Now we turn to the nucleon case. The nucleon scalar form factor G p (t) is defined by where M is the nucleon mass and u s (P ) is the Dirac spinor with spin s. From the Dirac structure, G p (t) is relevant to the nucleon helicity-flip amplitude at the high energy. Hence, we can apply the same methodology introduced in Sec. III E to investigate the large momentum transfer behavior of G p (t). It results in a factorization formula as follows: where the hard function has the following structure and H = H(y 1 ↔ y 3 ). Explicitly, we have where H Φ = H Ψ (1 ↔ 3) and the functions T i are defined by The above calculation reveals that the nucleon scalar form factor has the same power behavior as C g,q , i.e., G p (t) ∼ 1/(−t) 2 and hence P ↑ |F 2 |P ↓ ∼ 1/(−t) 3/2 . The convolution integral also suffers from the end point singularities, and the logarithm terms are expected.
With the above factorizations formula, we apply the methods introduced in Sec. IV to make numeric estimates on the gluon trace anomaly transition matrix for the nucleon, P ↑ | β(g) 2g F 2 |P ↓ . Similarly, the soft scale Λ c is introduced to regulate the end point singularities that appear in the hard part. In Fig. 8, we present the result in the range 0.35 GeV 2 < |t| < 5GeV 2 with three models of nucleon DAs in Table I P ↑ | β(g) 2g F 2 |P ↓ with respect to −t within different models of nucleon DAs.
negative in the large−|t| region. On the other hand, it is known that at zero momentum transfer, it yields the proton quantum anomalous energy and is positive [12]. That means that the transition matrix element should change sign when |t| varies from the low−|t| region to the high−|t| region. By inspecting the results obtained from the explicit models of nucleon DAs, one can see the sign change comes from the end point contributions, for example, in the asymptotic-DA model, where the constant term and the double-logarithm term in the square bracket have negative contributions to the transition matrix, while the single logarithm has positive contributions and dominate in the low-|t| region. Although the numerical results above are based on the perturbative analysis and in principle, only make sense in the large momentum transfer region(|t| Λ 2 QCD ), we expect the soft end-point contributions should be one of the origins of the sign change.

VI. CONCLUSION
In this paper, we have performed a detailed perturbative analysis of the gravitation form factors for nucleon at large momentum transfer, which results in the factorizations formulas of the GFFs in terms of the twist-3 or twist-4 nucleon distribution amplitude. We derived the hard coefficients at the leading order explicitly.
For A q,g form factors, we have demonstrated that they are related to the helicity-conserved amplitude and the quark and gluon contributions both scale 1/(−t) 2 . We have also computed B q,g , C q,g , C q,g form factors separately from the helicity-flip amplitude. Multiple consistent checks have been performed. We found that C q,g , B g,q are power suppressed as compared to the A p q,g . C q,g scale as A p q,g due to parameterization and the total C vanishes as expected. Because of the end point singularities, the B q,g , C q,g , C q,g form factors receive additional logarithmic contributions.
We have also computed the scalar form factors P | β(g) 2g F 2 |P for pion and proton and found out that both are negative at large (−t). It is understood that they are related to the hadron masses at zero momentum transfer. Therefore, there must be a sign change at lower (−t). This may lead to a nontrivial mass distributions of hadrons.
Based on the models for the twist-3 and twist-4 distribution amplitudes, we have performed numeric estimates for the quark and gluon GFFs in the perturbative region of large momentum transfer. It will be interesting to test these predictions in future experiments.
where the light-cone coordinates a ± = (a 0 ± a 3 )/ √ 2 are used. Suppose the nucleon state is moving along z direction with the momentum (P + , P − , 0 ⊥ ), the twist-3 nucleon distribution amplitude Φ 3 is defined by The twist-4 nucleon distribution amplitudes Φ 4 and Ψ 4 are defined by Hereû,d represent the u, d quark fields respectively and the gauge links that connect the fields are implied. C = iγ 2 γ 0 is the charge conjugation matrix. a, b, c refer to the color of the quark fields. M is the nucleon mass. The label [dx] denotes the measure dx 1 dx 2 dx 3 δ(x 1 + x 2 + x 3 − 1). In the Dirac representation, the Dirac spinors with definite helicity are taken as u ↑ (p) = 1 in which P and P are the initial and final nucleon momenta respectively with P 2 = P 2 = M 2 , s µ denote the covariant spin vector of the nucleon, and ψ q is the quark field of flavor q. L is the gauge link in the fundamental representation: where t a is the SU (3) generator in the fundamental representation. In the definition of GPD,P = (P + P )/2 is the average momentum, ∆ = P − P is the momentum transfer with t = ∆ 2 . The skewness parameter ξ is defined as ξ = P + −P + P + +P + , and the light-cone coordinates a ± = (a 0 ± a 3 )/ √ 2 are used.
Similar to the GFFs, the quark GPDs can also be factorized at the large −t limit. For nucleon H q , it is given by the nucleon helicity-conserved amplitude and its factorization has been reported in [92]. With the relation in Eq.(B1), we have checked our result on the quark A-GFF is consistent with theirs. On the other hand, the quark GPD E q is related to the nucleon helicity-flip amplitude, and its behaviors at large −t was unknown in the literature.
Following the discussions in the Sec.(III E), we derive the factorization theorem for the nucleon GPD E q at large momentum transfer. The formula reads where the hard perturbative coefficients E Ψq and E Φq have the following form: with u, d denote the quark flavors. At leading order perturbation theory, the explicit calculation yields where T iΨ and T iΦ are given by