Hadron tomography by generalized distribution amplitudes in pion-pair production process $\gamma^* \gamma \rightarrow \pi^0 \pi^0 $ and gravitational form factors for pion

Hadron tomography can be investigated by three-dimensional structure functions such as generalized parton distributions (GPDs), transverse-momentum-dependent parton distributions, and generalized distribution amplitudes (GDAs). Here, we extract the GDAs, which are $s$-$t$ crossed quantities of the GPDs, from cross-section measurements of hadron-pair production process $\gamma^* \gamma \rightarrow \pi^0 \pi^0$ at KEKB. This work is the first attempt to obtain the GDAs from the actual experimental data. The GDAs are expressed by a number of parameters and they are determined from the data of $\gamma^* \gamma \rightarrow \pi^0 \pi^0$ by including intermediate scalar- and tensor-meson contributions to the cross section. Our results indicate that the dependence of parton-momentum fraction $z$ in the GDAs is close to the asymptotic one. The timelike gravitational form factors $\Theta_1$ and $\Theta_2$ are obtained from our GDAs, and they are converted to the spacelike ones by the dispersion relation.From the spacelike $\Theta_1$ and $\Theta_2$, gravitational densities of the pion are calculated. Then, we obtained the mass (energy) radius and the mechanical (pressure and shear force) radius from $\Theta_2$ and $\Theta_1$, respectively. They are calculated as $\sqrt {\langle r^2 \rangle _{\text{mass}}} =0.32 \sim 0.39$ fm, whereas the mechanical radius is larger $\sqrt {\langle r^2 \rangle _{\text{mech}}} =0.82 \sim 0.88$ fm. This is the first report on the gravitational radius of a hadron from actual experimental measurements. It is interesting to find the possibility that the gravitational mass and mechanical radii could be different from the experimental charge radius $\sqrt {\langle r^2 \rangle _{\text{charge}}} =0.672 \pm 0.008$ fm for the charged pion.

Hadron tomography can be investigated by three-dimensional structure functions such as generalized parton distributions (GPDs), transverse-momentum-dependent parton distributions, and generalized distribution amplitudes (GDAs). Here, we extract the GDAs, which are s-t crossed quantities of the GPDs, from cross-section measurements of hadron-pair production process γ * γ → π 0 π 0 at KEKB. This work is the first attempt to obtain the GDAs from the actual experimental data. The GDAs are expressed by a number of parameters and they are determined from the data of γ * γ → π 0 π 0 by including intermediate scalar-and tensor-meson contributions to the cross section.
Our results indicate that the dependence of parton-momentum fraction z in the GDAs is close to the asymptotic one. The timelike gravitational form factors Θ1 and Θ2 are obtained from our GDAs, and they are converted to the spacelike ones by the dispersion relation. From the spacelike Θ1 and Θ2, gravitational densities of the pion are calculated. Then, we obtained the mass (energy) radius and the mechanical (pressure and shear force) radius from Θ2 and Θ1, respectively. They are calculated as r 2 mass = 0.56 ∼ 0.69 fm, whereas the mechanical radius is larger r 2 mech = 1.45 ∼ 1.56 fm. This is the first report on the gravitational radius of a hadron from actual experimental measurements. It is interesting to find the possibility that the gravitational mass and mechanical radii could be different from the experimental charge radius r 2 charge = 0.672 ± 0.008 fm for the charged pion. For drawing a clear conclusion on the GDAs of hadrons, accurate experimental data are needed, and it should be possible, for example, by future measurements of super-KEKB and international linear collider. Accurate measurements will not only provide important information on hadron tomography but also possibly shed light on gravitational physics in the quark and gluon level.

I. INTRODUCTION
Internal structure of hadrons has been investigated in terms of form factors and parton distribution functions (PDFs). Now, the field of hadron tomography, namely hadron-structure studies by three-dimensional (3D) structure functions [1][2][3][4], is one of fast developing areas in particle and nuclear physics. The 3D structure functions contain information on both the form factors and the PDFs, and they are ultimate quantities for understanding the nature of hadrons from low to high energies. Furthermore, it is essential to investigate the 3D structure of the nucleon for understanding the origin of nucleon spin because orbital angular momenta of partons could play an important role. The 3D structure functions could be also useful for clarifying internal quark-gluon configurations of exotic-hadron candidates [4].
Among the 3D structure functions, generalized parton distributions (GPDs) [1,2] and transverse-momentumdependent parton distributions (TMDs) [3] have been investigated extensively in recent years. We now have crude idea on these distributions. There are also generalized distribution amplitudes (GDAs) [1,4] as one of the 3D structure functions, and it is rather an unexplored field in comparison with the GPD and TMD studies. The GDAs can be obtained theoretically by the s-t crossing of the GPDs. Here, s and t are Mandelstam variables. Therefore, the GDA studies should also be valuable for the GPD understanding. In particular, both GPDs and GDAs can be expressed by common double distributions (DDs) with different Radon transforms as discussed later in Sec. II E. Therefore, the GDA studies are valuable also for understanding the GPDs through the DDs and simply by the s-t crossing.
The GDAs are key quantities for probing 3D structure of hadrons by timelike processes. In addition, one of the other important advantages of the GDAs is that 3D tomography is possible in principle for exotic-hadron candidates [4] because they can be produced in a pair in the final state, whereas no stable exotic hadron exists as a fixed target for measuring their GPDs and TMDs. The constituent counting rule can be used for identifying the number of elementary constituents in exotic hadron candidates at high energies. We should be able to distinguish exotic multiquark states from the ordinary qq and qqq ones by the counting rule [4,5]. Furthermore, form factors contained in the GDAs should provide information whether exotic hadron candidates are diffuse molecular states or compact multiquark ones [4].
Another advantage is that the GDAs and GPDs contain information on form factors of the energymomentum tensor so that the gravitational-interaction radius can be investigated. Although the root-meansquare charge radii are well known for the nucleons, the gravitational radius has never been measured experimentally. We try to extract the gravitational-interaction sizes, namely mass and mechanical radii, from existing experimental data in this work. Of course, the gravitational interactions are too weak to be measured directly for hadrons and elementary particles, such as quarks and gluons, "usually" in accelerator experiments, and there is no reliable quantum theory for the gravitational interactions at this stage. Nonetheless, it is interesting that the hadron tomography studies can access the gravitational information in hadrons through the energy-momentum tensor.
Fortunately, the Belle collaboration recently reported the cross sections for the pion-pair production in twophoton process γ * γ → π 0 π 0 at KEKB with various kinematical conditions [6,7]. It is our research purpose of this paper to extract the pion GDAs from the Belle measurements. Our studies should be the first attempt to extract any hadron GDAs from actual experimental measurements. Now, other hadron production processes γ * γ → hh are being analyzed in the Belle collaboration, so that other GDAs can be extracted in future. Furthermore, the KEKB accelerator has just upgraded and accurate measurements are expected in future for the two-photon processes. The two-photon processes have been used for investigating existence and properties of new hadrons in electron-positron annihilation reactions [8]. The same two-photon processes should be possible at the future international linear collider [9], and the GDAs will be investigated in the PANDA project [10]. This work is merely the first step for determining the GDAs; however, much progress is expected in the near future.
In this article, the generalized TMD (GTMD) or the Wigner distribution is explained first as a generating function for the 3D structure functions in Sec. II. Then, the GPDs and GDAs are introduced, and the form factors of energy-momentum tensor in the GDAs are explained in connection with the gravitational radii. Next, our theoretical formalism is developed for the γ * γ → π 0 π 0 cross section and the pion GDAs in Sec. III. The cross section of γ * γ → π 0 π 0 is expressed in terms of the GDAs. For extracting the GDAs from the experimental data, we introduce a parametrization of the GDAs, which are then determined by the analysis of the Belle measurement. Our analysis method is described in Sec. III, results are shown in Sec. IV. Finally, our studies are summarized in Sec. V.

II. THREE-DIMENSIONAL STRUCTURE FUNCTIONS OF HADRONS
The 3D structure of hadrons becomes one of hot topics in hadron physics, and it can be investigated by the GPDs, TMDs, and GDAs. First, we explain the Wigner phase-space distribution and the GTMD in Sec. II A as generating functions for form factors, PDFs, and the 3D structure functions. Then, we discuss the details of the GPDs and GDAs which are relevant to our studies including their relations in Secs. II B, II C, and II D. Both GPDs and GDAs are expressed by double distributions through Radon transforms as explained in Sec. II E. The GDAs are related to the timelike form factors of the energymomentum tensor, and then the spacelike gravitational form factors and radii are explained in Sec. II F.

A. Wigner distribution and three-dimensional structure functions
The 3D structure functions originate from the generating function, called the Wigner distribution, which is a phase space distribution W ( r, k ) expressed by the space coordinate r and momentum k. In the classical limit of → 0, it becomes the δ function δ(H( r, k )−E), which is the classical trajectory in the phase space. Therefore, its delocalization indicates quantum effects, and the Wigner function contains full information for describing quantum systems.
For the nucleon, the Wigner distribution was originally defined in Ref. [11] as the 6-dimensional phase-space distribution W (x, k T , r ), where x is the Bjorken scaling variable and k T is the transverse momentum. However, it was defined in a special Lorentz frame, so that a new definition was proposed in the infinite momentum frame [12] to express it by 5-dimensional phase-space distribution W (x, k T , r T ). It is equal to the ∆ + = 0 (ξ = 0) limit of the generalized transverse-momentum-dependent parton distribution (GTMD) [13].
In Fig. 1, relations of the GTMD and the Wigner distribution to the form factor, PDF, and 3D structure functions are shown [11][12][13] by integrating the GTMD by various kinematical variables. The form factors and the PDFs have been investigated until recently, and now the nucleon-structure studies focus on 3D structure functions, the GPDs and TMDs. However, there are few recent research activities on the GDAs which have a close connection to the GPDs by the s-t crossing.

B. Generalized parton distributions
The GPDs have been investigated by deeply virtual Compton scattering (DVCS) process as shown in Fig. 2 and also by meson and lepton-pair production processes [14,15]. In addition, there are possibilities to study them at hadron-beam facilities, for example, by N + N → N + π + B where N and B are the nucleon and baryon [16] and by an exclusive Drell-Yan process π − +p → µ + µ − +n [15,17,18]. Here, we explain the definition of the GPDs by using the DVCS process (γ * + h → γ + h) because its s-t crossing is the two-photon process (γ * + γ → h +h) which is analyzed in this work in terms of the GDAs.
We define kinematical variables for expressing the GPDs of the nucleon. The initial and final momenta of the nucleon are p and p ′ , respectively, as shown in Fig. 2, and they are q and q ′ for the photon. Then, their average momenta and the momentum transfer are given as [1,2,4,19] Expressing the momentum squared quantities as Q 2 = −q 2 andQ 2 = −q 2 , we define the Bjorken scaling variable x, momentum-transfer-squared t, and the skewdness parameter ξ as If the kinematical condition Q 2 ≫ |t| is satisfied, the skewdness parameter is expressed by the lightconecoordinate expression as The lightcone notation is given by a = (a + , a − , a ⊥ ) with a ± = (a 0 ± a 3 )/ √ 2 and the transverse vector a ⊥ . Then, the momenta are expressed as by using the relation (p + ) 2 , Q 2 ≫ M 2 , |t|. The scaling variable x is the lightcone momentum fraction carried by a quark in the nucleon, whereas the skewdness parameter ξ or the momentum ∆ indicates the momentum transfer from the initial nucleon to the final one or the momentum transfer between the initial and final quarks. The cross section of the DVCS γ * h → γh can be factorized into the hard part of quark interactions and the soft one expressed by the GPDs as shown in Fig. 2 if the kinematical condition is satisfied. Here, Λ QCD is the QCD scale parameter. The GPDs for the nucleon are defined by off-forward matrix elements of quark and gluon operators with a lightcone separation, and quark GPDs are defined by Here, q(y/2) is the quark field, M is the nucleon mass, and σ αβ is given by σ αβ = (i/2)[γ α , γ β ]. The functions H q (x, ξ, t) and E q (x, ξ, t) are the unpolarized GPDs of the nucleon, and there are also gluon GPDs H g (x, ξ, t) and E g (x, ξ, t) defined in a similar way [2]. To be precise, the link operator needs to be introduced in the left-hand side of Eq. (6) to satisfy the color gauge invariance. In this article, it is simply ignored. The advantages of the GPDs are that they contain both longitudinal momentum distributions for partons and transverse form factors. In fact, the GPDs H q (x, ξ, t) become unpolarized PDFs for the nucleon in the forward limit (∆, ξ, t → 0): where θ(x) is the step function, θ(x) = 1 for x > 0 and θ(x) = 0 for x < 0. Their first moments become Dirac and Pauli form factors F 1 (t) and F 2 (t), respectively: Another important feature, actually the most important for high-energy spin physicists, of the GPDs is that a second moment indicates a quark orbital-angularmomentum contribution (L q ) to the nucleon spin: because we know the quark contribution ∆q + = ∆q + ∆q from polarized charged-lepton DIS measurements. The GPDs have been mainly investigated for the nucleon. However, since the pion GDAs are investigated in this work and they are related to the pion GPDs by the s-t crossing, we also show the definition of the pion GPDs in the same way with Eq. (6) for the nucleon [20]: The pion is a scalar particle, so that the function E q (x, ξ, t) does not exist.
In comparison with PDF parametrizations, such studies are still premature for the GPDs due to the lack of experimental information. The simplest idea is to use the factorized form into the longitudinal PDF q(x) and the transverse form factor F T (t, x) at x [21]. For example, it is expressed as at ξ = 0 for x > 0. Namely, the GPDs contain information on both the PDFs and the form factors as already shown by the sum rules in Eqs. (7) and (8).

C. Generalized distribution amplitudes
If we exchange the s and t channels in the Compton scattering in Fig. 2, it becomes the two-photon process γ * + γ → h +h in Fig. 3. The GDAs describe the production of the hadron pair hh from a qq or gluon pair. We explain kinematical variables for describing the twophoton process and the GDAs [1,[22][23][24][25][26] as shown in Fig. 3. The initial photon momenta are denoted as q and q ′ , the final hadron momenta are p and p ′ , P is their total momentum P = p + p ′ , and k and k ′ are quark and antiquark momenta. One of the photon is taken as a real one with q ′ 2 = 0, and another one should satisfy the condition so that the two-photon process is factorized into a hard part and a soft one in terms of the GDAs as shown in Fig. 3 [27]. Here, W 2 is one of the variables in the GDAs, and it is the invariant-mass squared W 2 of the finalhadron pair. It is also equal to the center-of-mass (c.m.) energy squared s: The second variable ζ indicates the lightcone momentum fraction for one of the final hadrons in the total momentum P as shown in Fig. 3: Here, θ is the scattering angle in the c.m. frame of the final hadrons with the momentum assignments: and β is the hadron velocity defined by with the final-hadron mass m h . The third variable z is the lightcone momentum fraction for a quark in the total hadron-pair momentum P , and it is defined by The GDAs are expressed by these three variables, z, ζ, and W 2 = s. The quark GDAs are defined by the matrix element of the same operators used in defining the GPDs in Eq. (6) between the the vacuum and the hadron pair: We use the notation Φ hh q for one specific quark (q) without the summation over the quark flavor. Here, the kinematical range of z is 0 ≤ z ≤ 1, whereas the variable z ′ = 2z − 1 is often used with the same notation z (or x) in the range −1 ≤ z ′ ≤ 1 for the distribution amplitude as explained in Ref. [18]. However, because many articles of the GDAs use the notation z in the range 0 ≤ z ≤ 1, we follow this convention in this work. The expression e i(2z−1) P + y − /2 h(p)h(p ′ ) | ψ(−y/2)γ + ψ(y/2) | 0 is sometimes written by the equivalent one as e −iz P + y − h(p)h(p ′ ) | ψ(y)γ + ψ(0) | 0 . Furthermore, the gauge link should be introduced in the nonlocal operator to satisfy the color gauge invariance; however, it is simply neglected in this paper. There are sum rules for the quark GDAs of the isospin I = 0 two-meson final states [23,24]: where M h 2(q) is the momentum fraction carried by flavor-q quarks and antiquarks in the hadron h (note: total quark fraction q M h 2(q) ). As shown in Eq. (44), this integral is expressed by the energy-momentum tensor of a quark, so that the right-hand-side of Eq. (19) should be described by the form factors of the energy-momentum tensor at finite W 2 [28]. There are recent theoretical studies on the energy-momentum tensor for the nucleon [29] and on its lattice QCD estimate [30]. In general, there are two energy-momentum tensor form factors for the pion [31,32], and they are explained in Secs.II F and III H.
Since the GDAs contain intermediate-meson contributions as explained in Sec. III E, the second sum of Eq. (19) should be a complex value at finite W 2 . There are resonance terms and the continuum one which contains a quark part of the form factor F h q (W 2 ) defined in Eq. (94). The explicit expression is shown later in Eqs. (124) and (125) for analyzing actual experimental data. Therefore, our studies can suggest the optimum form factor F h q (W 2 ) of the energy-momentum tensor for the continuum part of the hadron h, and they are related to the size of gravitational interaction. The gravitational radii of a hadron are discussed in more details in Sec. II F. The sum rule of Eq. (19) was derived for the kinematical point of W 2 = 0 [23,24], and then it was considered even at finite W 2 as the form of form factor of the energy-momentum tensor [28]. However, since there are two gravitational form factors for the pion in general, a relation between the GDAs and the form factors is newly derived in Sec. III H of this article.
The GDAs are defined for the hadron-antihadron system, so that they satisfy the charge-conjugation invari- ance [2]: where C is the charge-conjugation operator. We may note that the gluon GDA should satisfy the condition due to the translational invariance in defining the gluon GDA and C invariance. As shown in Fig. 4, the gluon GDA contributes to the two-photon cross section as a next-to-leading order term [24,25], so that it is neglected in our current leading-order analysis.

Generalized distribution amplitudes for pions
The pion GDAs are investigated in this work, and there are two notation types for them. One is the representation based on the C-parity eigenstates, and the other is by the isospin. In order to avoid confusion, we explain them here in details.
First, we consider possible two-pion states. Denoting I for the isospin, we have the I = 1 ππ state which is antisymmetric under the exchange of the pions. On the other hand, the I = 0 and I = 2 ππ states are symmetric: Since the C parity of γ * γ is even, the ππ state needs to satisfy C = (−1) L+S = (−1) L = +1 with S = 0. Therefore, L should be even. The Pauli principle indicates so that the isospin states should be I = 0 or 2. The GDAs are defined in Eq. (18) by the the matrix element of the vector-type nonlocal operator. Since the isospin ofqq is 0 or 1, the only possible choice for the ππ isospin is I = 0. In this way, only the ππ sates allowed in the γ * γ process should have I = 0 with L = even numbers (0, 2, · · · ). For the actual pion GDAs which are investigated in this work, we may express them as C-parity eigenstates [24] Φ ππ(±) q (z, ζ, where (±) indicates the C parity. Therefore, the π + π − GDAs are given by , and the C-even part satisfies Φ (1 − z, ζ, W 2 ). The π 0 π 0 GDAs contain only the C-even function: Then, the isospin invariance leads to the relations between the u-and d-quark GDAs as Φ On the other hand, the isospin decomposition of the pion GDAs is discussed in Ref. [23] first by defining them as the twist-2 chiral-even amplitudes by where ψ is the quark field with u and d quark components for the isosinglet (isovector) GDA. The notationÎ is the identity matrix. They are expressed by the isoscalar and isovector GDAs as They satisfy the symmetry relations due to the charge conjugation: If the isospin-symmetry relations are satisfied for the pion GDAs, the isosinglet and isovector GDAs are related to the C even and odd GDAs as In this work of the π 0 π 0 production process, only the following isoscalar or C-even GDAs are involved in the cross section γ * γ → π 0 π 0 : where q indicates u or d. This function is parametrized and used for the analysis of π 0 π 0 production data later by using Eq. (65).

D. Relation between GPDs and GDAs
As obvious from the diagrams of the DVCS and twophoton process in Figs. 2 and 3, respectively, the GPDs and GDAs are related with each other by the s-t crossing as long as the factorization conditions are satisfied. Namely, the scale Q 2 should be large enough for the factorization: By the s-t crossing, the final hadronh with the momentum p ′ becomes the initial hadron h with p, which indicates the momentum changes from p and p ′ in the GDAs to p ′ and −p in the GPDs. It means that both variables are related by the relations [1,24].
and then the GDAs are GPDs are related to each other by The physical regions of the kinematical variables are However, the relation of Eq. (32) indicates that the physical GDAs do not necessary correspond to the physical regions in Eq. (33) of the GPDs: Namely, the GDAs could lead to the unphysical kinematical regions, |x| > 1, |ξ| > 1, and t > 0, of the GPDs. Equation (32) also indicates the relation |ξ| ≥ |x|, which is called as the Efremov-Radyushkin-Brodsky-Lepage (ERBL) region. The ERBL region of the GPDs can be investigated, for example, by the hadronic reaction N + N → N + π + B [16]. However, GDA studies will provide another information on the ERBL GPDs although it is in the unphysical region of t > 0.

E. Radon transforms for GPDs and GDAs by using double distributions
We explained definitions and basic properties of the GPDs and GDAs. They are related with each other by the s-t crossing. The studies of the GDAs should be valuable for the GPD studies and vice versa. In fact, both GPDs and GDAs are expressed by the common double distributions (DDs) by different Radon transforms. The Radon transform is defined in n dimensions for an arbitrary function f (x) by [33] where x is the n-dimensional space coordinate [x = (x 1 , x 2 , · · · , x n )] and ξ is the unit vector in n dimensions [ξ = (ξ 1 , ξ 2 , · · · , ξ n )]. Because of the δ function, the integral is over the n − 1-dimensional plane constrained by p = ξ · x.
Using this Radon transform, we can express the GPDs and GDAs in terms of double distributions (DDs), F q (β, α, t) and G q (β, α, t), defined by the matrix element [1,20] for the scalar hadron h like the pion. The kinematical support region is given by |β| + |α| ≤ 1 for the DDs.
Using the Radon transform, we can express the GPDs in terms of these DDs as Namely, the GPDs are obtained by integrating the DDs over the slight line x = β + ξα as shown in Fig. 5. The parton distribution functions (PDFs) are obtained as a special case of this integral over the vertical line in Fig. 5 with the constraint of the forward limit (t = 0), and they are expressed as There are similar relations of the gluon DD to the gluon GPDs and PDF [1].
As just an example, we introduce a simple parametrization for the DDs F q (β, α), which are expressed by the corresponding PDF q(β) multiplied by a profile function h q (β, α) as [34] Kinematical support region of the double distributions and integral paths for obtaining the GPDs, GDAs, and PDFs.
The profile function may be expressed as The matrix element associated with the GDAs is also expressed in the same way by the DDs as [1,20] Then, the GDAs can be expressed by the DDs as which indicates that the GDAs are obtained by the Radon transform along the different line 1 − 2z − (1 − 2ζ)β + α = 0 as shown in Fig. 5. We found that both GPDs and GDAs can be expressed by the DDs. Therefore, experimental measurements of the GDAs should be valuable also for the GPD studies through the determination of the DDs and vice versa. In particular, the GDAs correspond to specific kinematical regions of the GPDs as explained in Sec. II D. These investigations from the direction of the GDAs could be supplementary to the direct GPD studies. Furthermore, it is the advantage of the GDAs that exotic hadron GDAs can be measured in future, whereas their GPDs cannot be studied experimentally because there is no stable exotic-hadron target. Considering these merits, we believe that our GDA project should be important for future developments on hadron tomography not only for ordinary hadrons such as the nucleons and pions but also for exotic-hadron candidates.

F. Timelike form factors of energy-momentum tensor and gravitational-interaction radii
The GPDs and GDAs are measured in the DVCS and two-photon processes which are, of course, electromagnetic interaction processes. However, their studies could also probe an aspect of gravitational interactions with quarks and gluons. In order to understand this fact, we explain it by taking the quark GPD and GDA definitions. As given in Eqs. (6) and (18), the GPDs and GDAs are defined by the same non-local vector operator. For the GDAs, its moments multiplied by the momentum factor 2(P + /2) n are expressed by the derivatives as [1] For n = 2, this operator is the energymomentum tensor of a quark, and it is a source of gravity, whereas it is the vector-type electromagnetic current for n = 1.
As shown in Fig. 6, (a) the electromagnetic interaction is described by the vector currentqγ µ q, (b) the weak interaction is characterized by the vector minus axialvector current γ µ (1 − γ 5 ), and (c) the gravitational one is by the tensor interaction given byqγ µ ∂ ν q for a quark. In Eq. (43), the GPDs and GDAs contain this factor as the energy-momentum tensor of a quark. The charge radius of the proton is measured by elastic electron scattering in the form of the electric form factor through the photon exchange process (a). In the similar way, the gravitational radius should be measured by the graviton exchange process (c) in principle. However, it is impossible to do an actual scattering experiment directly at accelerator facilities for the gravitational interaction due to the ultra-weak interaction nature. On the other hand, it is possible to access such physics through the GPDs and GDAs. Therefore, the gravitational radii of hadrons are measurable quantities, although it may be somewhat surprising that a different physics aspect can be investigated through the electromagnetic processes.
In this way, the following integral of the quark GDAs over the variable z is related to a matrix element of the quark energy-momentum tensors T µν q [23,24,28]: and there is a similar equation on the gluon matrix element. The quark energy-momentum tensor is given by where D µ is the covariant derivative D µ = ∂ µ − igλ a A a,µ /2 defined by the QCD coupling constant g and the SU(3) Gell-Mann matrix λ a . The notation X (µν) is given by the symmetric combination X (µν) ≡ (X µν + X νµ )/2. Here, T ++ q indicates the lightcone ++ components as expressed in Eq. (43), so that it is specifically given by T ++ q = (T 00 q + T 03 q + T 30 q + T 33 q )/2. These equations indicate that the GDAs probe the energymomentum tensors of quarks and gluons, in the same way as the GPDs, in the timelike process.
In an isolated system, the energy-momentum tensor is conserved ∂ µ T µν = 0. However, if there is an external force and gravity, it satisfies [35] 1 where G ν is the energy-force density, g is defined by the metric tensor g µν as g = det(g µν ), and Γ ν µλ is the affine connection tensor. The second term on the right-hand side is the gravitational-force density. As the electromagnetic interaction and weak interaction are characterized theoretically by the vector and vector minus axialvector operators,qγ µ q andqγ µ (1 − γ 5 )q, respectively for quarks, the gravitational interaction is characterized by the energy-momentum tensor T µν . Namely, the energymomentum tensors of quarks and gluons are sources of gravitational interactions. Now, the 3D structurefunction studies are getting popular in hadron physics, and this tensor appears in the 3D structure functions, in particularly in the GPDs and GDAs as illustrated in Fig. 6. For example, the GDAs probe the 3D structure of a hadron in the form of the timelike form factors. The GDAs are related to the energy-momentum tensor in Eq. (44), so that they probe the gravitational interaction, for example, as the form factors of energy momentum tensor. These form factors are explicitly defined later in Eqs. (57) and (116).
The GPDs and GDAs contain information on spacelike and timelike form factors, respectively. For example, the simple parametrization of the GPDs is given in Eq. (11) expressed as the longitudinal PDF multiplied by the twodimensional transverse form factor. In general, the twodimensional transverse charge density ρ h T (r ⊥ ) is given by the Fourier transform of the spacelike electric form factor of a hadron h as where J 0 is the Bessel function. The two-dimensional transverse root-mean-square (rms) radius is then given by . (48) The transverse form factors of the energy-momentum tensor are calculated by using a simple parametrization for the GPDs of the proton, and the results indicate that they could be different from charge form factor [36].
In the three-dimensional case, the charge density and the form factor are related with each other by where j 0 is the spherical Bessel function. The rms radius is obtained by For timelike form factors probed by the e + e − or γ * γ reactions, we can relate them to the spacial distributions by using the dispersion relation. Considering that singularities of the form factor F h (t) is in the positive real t axis from 4m 2 h , we can express the t-channel form factor by the dispersion integral over the real positive t (≡ s) as [37,38]: Namely, the t-channel form factor F h (t) can be calculated from the s-channel one F h (s). Then, using Eqs. (47) and (51) together with consideration on the constituent-counting rule in the asymptotic region [4], we have [38] where K 0 is the modified Bessel function of the second kind. However, the imaginary part of the form factor, namely its phase, is not available from the measurement of γ * → hh because its cross section is proportional to |F h (t)| 2 and a theoretically model-dependent input is needed for estimating the spacial charge distribution from the measurement on the timelike form factor. In Ref. [38], the Gounaris-Sakurai amplitude [39] is used for ImF π (t) to obtain the transverse charge radius r 2 ⊥ π ch = 0.53 fm, which corresponds to the threedimensional one r 2 π ch = 1.5 r 2 ⊥ π ch = 0.42 fm 2 . Here, "ch" indicates the electric charge. This value is comparable to the πe scattering measurement value r 2 π ch = 0.439 ± 0.008 fm 2 [40] for the charged pion. These results are for electric charge radii probed by electromagnetic interactions, whereas we investigate gravitational radii for hadrons, particularly the pion in this work, by using the GDAs in the two-photon process γ * γ → hh.
The three-dimensional density is calculated by using Eqs. (49) and (51) as The three-dimensional rms radius is also obtained by using Eqs. (50) and (51) as The spacelike gravitational form factors Θ 1 (t) and Θ 2 (t) are defined by the energy-momentum tensor T µν [28,31,32]. In the GPD and GDA studies [28], other notations A(t) and B(t) are often used. Here, A and B are used for expressing other quantities, so that we use the notations Θ 1 (t) and Θ 2 (t) for the gravitational form factors. In the spacelike process, they are defined by for a quark q. Here, the momenta are defined by P = p + p ′ and q = p ′ − p. We defined the form factors and the energy-momentum tensor for one quark type (namely, flavor-q quark and antiquark) in order to avoid confusions. In Ref. [28], the form factors are expressed by A and B, and they are related to Θ 1 (t) and Θ 2 (t) by As discussed above Eq. (31), the variables (p, p ′ ) (GPD) in the t channel is changed for (−p ′ , p) (GDA) in the s channel by the s-t crossing. Then, using the momentum notations P = p + p ′ and ∆ = p ′ − p, we obtain the definition of the timelike form factors from Eq. (55) as From Eq. (44) and this definition, we can evaluate the gravitational form factors for the pion.

III. THEORETICAL FORMALISM
We explain the cross section for the two-photon process γ * γ → π 0 π 0 to express it in terms of the GDAs in Sec. III A. First, the situation of the pion distribution amplitude (DA), instead of the GDAs, is explained in Sec. III B, and Q 2 evolution of the DA and the GDAs are discussed in Sec. III C. The ζ dependence of the GDAs is introduced in Sec. III D. Then, the parametrization of the GDAs is introduced in Sec. III E to determine them from experimental data. Contributions from f 0 and f 2 resonances are included in the analysis, and coupling constants for the resonances are explained in Sec. III F, and the Q 2 scale dependence of such resonance contributions is discussed in Sec. III G. In Sec. III H, the relations between the gravitational form factors and the GDAs are derived.
A. Cross section for the two-photon process The pion-pair production process γ * γ → π 0 π 0 is shown in Fig. 7, and its cross section is written by the matrix element M as [41] where one of the initial photons is taken on mass shell (q ′ 2 = 0). The matrix element M(γ * γ → π 0 π 0 ) is given by the hadron tensor T µν as by the photon polarization vector ǫ µ (λ) and the electromagnetic current J em µ (y). In obtaining the total cross section, the cross section should be divided by two due to two identical particles in the final state to avoid the double counting. Alternatively, the cross section is integrated over the half solid angle, in stead of the factor 1/2, for calculating the total cross section. In any case, differential cross sections are discussed in this paper, so that such a factor is not needed.
We define the helicity amplitudes A ij by If the kinematical condition Q 2 ≫ W 2 , Λ 2 is satisfied, the two-photon process can be factorized into the hard part (γ * γ → qq) and the soft part (qq → π 0 π 0 ) as shown in Fig. 3. In the Breit frame, q is taken along the z axis. Introducing two timelike vectors n = (1, 0, 0, 1)/ √ 2 and n ′ = (1, 0, 0, −1)/ √ 2, we express the photon and quark momenta as At large Q 2 , the hadron tensor can be expressed by the factorized form as The first part describes the process γ * γ → qq of Fig. 3, and the second one does the soft process qq → π 0 π 0 . For the term q b (y)q a (0) in this equation, we use the Fierz identity where the first two terms and the last one are the leading twist terms, while the third and fourth ones are twist-3 terms. Since the trace of an odd number of γ λ is zero, only the first two terms survive. However, the second term is the axial-vector current, which cannot exist for π 0 π 0 state due to the party invariance. After all, only the first term contributes to the hadron tensor.
In the leading order of the running coupling constant α s , the gluon GDA contribution is neglected and the hadron tensor can be expressed by the quark GDAs by calculating the hard part of Eq. (61) as [4,23,24] where g µν T is defined by The hadron tensor T µν is generally written by the product of the two electromagnetic currents in Eq. (59). In the leading twist, it is expressed by the matrix element of the vector current as given by the GDAs Φ π 0 π 0 q in Eq. (26) [24]. The situation is the same as the one in the hadron tensor W µν in the charged-lepton deep inelastic scattering as expressed in the twist expansion [42].
Since only the non-vanishing terms are ε where the relation A −− = A ++ is used due to party conservation. The gluon GDA contributes to the cross section through the amplitudes A ++ = A −− and A +− = A −+ in the next-to-leading order, so that these terms are suppressed by the factor of α s . There are also contributions from higher-twist amplitudes A 0+ and A 0− , which decrease as at least 1/Q because of a helicity flip [24,25]. The γ * γ → π 0 π 0 cross section is expressed by the GDAs in Eq. (65). In order to determine the GDAs from experimental data, we need to express the GDAs by a number of parameters, which are then determined by a χ 2 analysis of the data on dσ/d(cos θ). There a number of studies on the pion distribution amplitudes; however, it is the first attempt for the GDAs in comparison with actual experimental data. Before discussing an appropriate functional form of the GDAs, we explain the distribution amplitude (DAs), which are related to the z-dependent part of the GDAs. For example, the pion distribution amplitude Φ π (z) is related to the GDAs by [23] In our analysis of γ * γ → π 0 π 0 , we obtain Φ ππ(+) .
The function Φ π (z, µ) is the leading-twist distribution expressed by the longitudinal momentum fraction z of a valence quark in the pion and the renormalization scale µ of the bilocal operator. The µ dependence is described by the ERBL evolution equations [44]. It is normalized as and f π is the pion decay constant defined by In the asymptotic limit of µ → ∞, the pion distribution amplitude becomes as it becomes obvious from the Q 2 -evolution solution of Eqs. (78) and (81). At finite µ, it is generally expressed by using the Gegenbauer polynomials as Φ π (z, µ) = 6 z (1 − z) ∞ n=0,2,4,··· a n (µ) C 3/2 n (2z − 1), (71) where only the even terms contribute because the DA should satisfy the condition Φ π (1 − z, µ) = Φ π (z, µ) under the exchange z ↔ 1 − z. It corresponds to the exchange of q andq in the pion, and the momentum distribution carried by a quark or antiquark should be same under this exchange because of positive C-parity of the axial current. The Gegenbauer polynomials are The current situation of the pion DA is explained in Ref. [18]. Since the Gegenbauer polynomials are rapidly oscillating functions at large n and the coefficients a n (µ) are small for large µ, the n ≥ 4 terms could be neglected at this stage. As for the second coefficient a 2 , there are theoretical estimates by lattice QCD [45] and QCD sum rules [46][47][48][49][50][51] One of the well known functions was proposed by Chernyak and Zhitnitsky (CZ) to take a 2 (µ ≃ 0.5 GeV) = 2/3 as suggested by the QCD sum rule [46]: which is very different from the asymptotic form because it has a minimum at z = 0.5. There are also recent theoretical suggestions on different a 2 values [45,[48][49][50][51][52] and also a 4 and a 6 values [53]. In principle, the different pion DAs can be tested by experiments. The Belle measurements on the γ → π form factor are close to the asymptotic DA form [6], whereas the BaBar data have a different tendency in the sense that it is consistent with a 2 (µ = 2 GeV) = 0.22 [52]. Further measurements are needed to distinguish various theoretical DAs. We comment on a slightly different convention from ours in defining the distribution amplitude because it may be sometimes confusing in using the decay constant f π or f π / √ 2. In the Diehl's article of 2003 [1], the π 0 distribution amplitude is defined for one quark flavor as π 0 (p) | q(y) α q(0) β | 0 instead of the left-hand side of Eq. (67). Therefore, the decay constant definition of Eq. (69) We should note that there is a factor of √ 2 in this expression. However, this √ 2 is absorbed into the definition of distribution amplitude in his formalism so that the decay constant f π stays the same as ours.

C. Scale evolution of distribution amplitudes and generalized distribution amplitudes
If the kinematical condition Q 2 ≫ W 2 , Λ 2 QCD is satisfied, the process γ * γ → π 0 π 0 is factorized into the hard part H q,g and the soft one S q,g as shown in Fig. 8. Here, the final state X is π 0 for the DAs or π 0 π 0 for the GDAs. The hard part is calculated in perturbative QCD and the soft one is expressed by the DAs or the GDAs. The Q 2 evolution equations of the DAs and GDAs are described by calculating the hard part in perturbative QCD. Since both reactions (γ * γ → π 0 and γ * γ → π 0 π 0 ) have the same hard processes, the DAs and GDAs follow the same evolution equations, and their z and scale-µ dependencies are represented by the functions Φ q (z, µ) and Φ g (z, µ) in the following discussions of this subsection [Φ q = Φ π for the DAs, Φ q = Φ ππ(+) q and Φ g = Φ ππ g for the GDAs]. In order to describe the Q 2 evolution of the DAs and GDAs, we introduce the auxiliary quark and gluon func- 8. Factorization of γ * γ → X, X = π 0 for DAs or π 0 π 0 for GDAs, into the hard part Hq,g and the soft one Sq,g. For the isovector π, the right-hand-side process with the two-gluon intermediate state does not exist for the single π 0 production (X = π 0 ).
tions f Q and f G defined by [1,24] We introduce the variable τ defined by for describing the evolution from µ 2 0 to µ 2 as usually used in expressing the DGLAP evolution equations for the PDFs [54]. Here, β 0 = 11 − 2n f /3 and α s is the running coupling constant. Then, the evolution equations are expressed as where the matrix V is the kernel calculated in perturbative QCD. The one-loop contributions to this kernel is shown in Fig. 9. The one-loop kernels have been obtained as [1,24] V whereū andz are defined byū = 1 − u andz = 1 − z, and C F , T F , and C A are given by C F = (N 2 c − 1)/(2N c ), T F = 1/2, and C A = N c with the number of colors N c = 3. These equations are called ERBL evolution equations.
The integro-differential equations can be solved in the same way with solving the DGLAP (Dokshitzer-Gribov-Lipatov-Altarelli-Parisi) evolution equations [54] by using the anomalous dimensions (γ QQ n , γ QG n , γ GQ n , γ GG n ) FIG. 9. Leading contributions to the hard part Hq,g.
obtained from the kernel matrix V . From these anomalous dimensions, we define Then, the solution is written in terms of the Gegenbauer polynomials C a n and the anomalous dimensions as where the coefficients are given by with the factor g ± n = (γ ± n − γ QQ n )/(3γ QG n /n). The summations of Eq. (78) are taken for odd n (n = 1, 3, · · · ) in the C =even case.
In the n =odd summation, all the anomalous dimensions are positive except for γ − 1 = 0, so that the only the A − 1 terms survive in the scaling limit of µ 2 → ∞. Using the Gegenbauer polynomial C a 1 (x) = 2ax, we have the C =even (isoscalar) GPDs as Therefore, the z-dependent functional forms are uniquely given for the GDAs. This fact should be taken into account for parametrizing the GDAs. For the C = odd (isovector) GDAs, the n summation of Eq. (78) is for the even numbers (n = 0, 2, · · · ). In the scaling limit, only the n = 0 term survives and the GDAs become by using C a 0 (x) = 1. The above isovector GDAs have z dependence z(1 − z) which is the same as the ρ-meson (pion [23]) isovector DA of Eq. (70) in the scaling limit.

D. ζ dependence of generalized distribution amplitudes
The Q 2 evolution of the GDAs are calculated in perturbative QCD as shown in the previous subsection, and the z dependence is given by the Gegenbauer polynomials. The GDAs also depend on other two variables ζ and W 2 . Here, we discuss the ζ dependence. As shown in Fig. 3 and Eq. (14), the variable ζ indicates the momentum fraction for a produced pion in the final state and it is expressed by the polar angle (θ) of the pion. Therefore, we may expand the coefficients A n and A ′ n in terms of orthogonal polynomials, which could be taken as the Legendre polynomials P l : where n is odd (l is even) for C = +, and n is even (l is odd) for C = −. Here, the factor 6 comes in the similar way to the normalization of the pion DA as shown in Eqs. (68) and (70), the flavor number n f appears because of the flavor summation in Eq. (78), and l is the angular momentum of the final pion pair. In addition, the same equation exists for A ′ n (ζ, W 2 ) in terms of B ′ nl (W 2 ). The C invariance relations of the GDAs are given in Eqs. (20) and (21), so that the odd-l terms do not contribute to the C = + GDAs.
From the scale-dependence relations of Eq. (79), the coefficients B nl should follow the same relations: and a similar equation for B ′ nl (W 2 , µ). In the scaling limit µ → ∞, only the lowest terms survive in A n (ζ, W 2 ) and A ′ n (ζ, W 2 ), and we obtain [1,24] where the Legendre polynomial P 2 (x) is given by P 2 (x) = (3x 2 − 1)/2. Since the Legendre polynomial term is given by P 2 (2ζ − 1) = 1 − 6ζ(1 − ζ), the sum rule of Eq. (19) is satisfied if the coefficients satisfy the relation B 10 (W 2 = 0) = −B 12 (W 2 = 0), which is considered in the parametrization in the next subsection. This is the basic functional forms for z and ζ dependencies in the scaling limit. Next, we explain our actual parametrization for the GDAs by following the essence of these basic functional forms.

E. Expression of generalized distribution amplitudes
With the basic knowledge of the pion DA and GDAs, we need to express the GDAs by a number of parame-ters. In particular, the z dependence is given by Eqs. (80) and (81) in the scaling limit. Considering these functional forms, we express the GDAs with a number of parameters. First, we neglect the higher-order α s effects and higher-twist effects, so that the gluon GDA does not appear. Since π 0 π 0 -production data are analyzed in this work, only the C = even GDAs contribute to the cross section. The C = even function of Eq. (80) is z(1 − z)(2z − 1). Since the C = even isoscalar GDAs have − sign under the change z → 1 − z as given in Eq. (28), the same parameter α is assigned for the powers of z and 1 − z: Φ π 0 π 0 q (z) ∼ z α (1 − z) α (2z − 1). The 2z − 1 factor comes from the lowest Gegenbauer polynomial C 3/2 1 (2z − 1), which survives in the scaling limit. However, the detailed z dependence is not determined from the current data, so that we decide to take the lowest Gegenbauer polynomial form of 2z −1, supplementing by the phenomenological parameter α which will later appear to be close to the asymptotical value 1.
We use the following function for explaining the γ * γ → π 0 π 0 data at a fixed Q 2 value: where N α is the overall constant determined by the sum rule (19) as with the beta function B(a, b). The quark-momentum fraction factor M π 2(q) and the W 2 -dependent form factor F π q (W 2 ) are included in the coefficients B nl (W 2 ). The ζ dependence can be re-expressed by the angle θ defined in Eq. (15) as where the invariant-mass dependent functions B nl (W 2 ) and B nl (W 2 ) are related with each other by In the limit of W 2 → 4m 2 π ≃ 0, they are given by [23,24] B 12 (0) = 10 9 M π 2(q) , where M π 2(q) is the momentum fraction carried by the qflavor quarks and antiquarks in the pion ( n f q M π 2(q) ≃ 0.5). This equation is obtained by considering the forward limit of the GPDs and then the s-t crossing to relate the GPDs and GDAs, so that it should be a modelindependent relation. Then, the relation between B 10 (0) and B 12 (0) is studied in a soft-pion theorem, and it was obtained as [23,24] B 10 (0) = −B 12 (0).
Then, the W 2 dependence of B 10 (W 2 ) and B 20 (W 2 ) was studied at small W 2 as a possible constraint on the functional form of W 2 within a instanton model of QCD [28]. The gluon GDA does not contribute to the cross section because the higher-order and higher-twist terms are neglected in our analysis. However, as discussed in Sec. III C, it affects the Q 2 evolution. It will be shown in Figs. 13 and 14 that current Belle data are not accurate enough to probe the scaling violation. The quark GDAs are provided at a fixed Q 2 scale which is taken as the average Q 2 value (16.6 GeV 2 ) of the Belle data. Then, the Q 2 evolution is not taken into account in our analysis within the Belle-data range (8.92 ≤ Q 2 ≤ 24.25 GeV 2 ). Therefore, the gluon GDA does not contribute in our analysis.
There are two terms, which correspond to the angular momenta, l = 0 and 2, of the pion pair. There are intermediate meson contributions to the cross section for γ * γ → π 0 π 0 , so that the invariant-mass dependent factors B nl have imaginary parts expressed by the phase shifts δ l (W ): Here, we use the ππ phase shifts by Bydzovsky, Kaminski, Nazari, and Surovtsev [55]. There is also another study on the phase shifts in Ref. [56]. The relation of Eq. (89) indicates that theB nl (W 2 = 0) factors are given byB There are two types of contributions to B nl (W ). One is the continuum and the other is from the intermediate resonances expressed bȳ where M R is the resonance mass, Γ R is its width, and c R is a constant. The W 2 dependence of the continuum part of the pion GDAs is given by the form factor, which could be parametrized as [4] F π q (W 2 ) = Here, Λ is the cutoff parameter, which indicates the pion size, n is the number of active constituents according to the constituent-counting rule in perturbative QCD [57], and it is normalized as F h(q) (4m 2 π ) = 1. It is the continuum part of the timelike forms factor of the energy-momentum tensor. Here, the pion size means the gravitational-interaction size instead of the usual charge radius in electromagnetic interactions as explained in Sec. II F. The high-energy behavior of the form factor is given by the factor n, which is supposed to be n = 2 for the pion [4].

F. Resonance terms and their coupling constants
The resonance contributions are illustrated as the intermediate states in Fig. 10. Above 1 GeV of the invariant mass W , the intermediate KK and ηη channels contribute to the process. However, their contributions may not be as large as the pion ones, and they are not explicitly considered in this work. As seen in Fig. 10 in principle. However, these meson effects have minor effects on the cross section, hence they are not included in our analysis. The GDAs are defined by the matrix element of the nonlocal vector operator from the vacuum to the π 0 π 0 state, and it is expressed by three steps for describing the process with the intermediate f 0 meson [58]. First, the f 0 meson is produced from the vacuum, it propagates, and then it decays into the pion pair: In Sec III B, the pion distribution amplitude is defined by ψγ µ γ 5 ψ instead of one quark flavor one qγ µ γ 5 q. The above f 0 distribution amplitude is related to the one defined by ψγ µ ψ = (uγ µ u + dγ µ d)/ √ 2 as where q = u or d. The final 2π-decay part is simply the coupling constant written as and the first f 0 production part is expressed by the distribution amplitude for f 0 as discussed in the later part of this subsection. Then, the f 0 contribution to B 10 (W ) is written as so that its absolute value is given bȳ where the factor 5/3 comes from the convention difference in defining the distribution amplitude, namely the overall factor could be 30 or 18. This difference becomes 30/18 = 5/3. In the same way, the f 2 contribution is given bȳ where the different factor 10 M 2 f2 /9 comes from the tensor nature of f 2 in defining the coupling constant to 2π and also the decay constant [58,59]. In Ref. [58], the β 2 factor is included in Eq. (A26) of this paper.
As for the resonance terms, we use the W dependence of | B nl (W 2 )| in Eqs. (101) and (102) [59], although the resonance properties are also obtained by the Belle collaboration for the resonances f 0 (980) and f 2 (1270) [7]. In Refs. [59][60][61], the constants are f f2 = 0.101 GeV at Q 2 = 1 GeV 2 , M f2 = 1.275 GeV, and Γ f2 = 0.185 GeV for f 2 (1270), and the decay constant g f2ππ is defined by Here, the factor of 2 in 2/3 comes from the identical particles of two π 0 's, and the factor 1/3 does  Table III, the decay constants are listed at Q 2 = 16.6 GeV 2 , which is the average scale of the used Belle data. The value f f 2 = 0.101 at Q 2 = 1 GeV 2 corresponds to f f 2 = 0.0754 at Q 2 = 16.6 GeV 2 . from Γ(f 2 → π 0 π 0 ) = 1/3Γ(f 2 → ππ). As it will become clear in comparison with the actual measurements, the Belle data indicate a clear peak of f 2 (1270).
For the S-wave resonances of f 0 (500) and f 0 (980), we have M f0(500) = 0.475 GeV and Γ f0(500) = 0.55 GeV [62], the decay constant g f0ππ is defined by g f0ππ = (2/3)16πΓ(f 0 → ππ)M f0 with Γ(f 0 → ππ) = Γ f0 for both f 0 (500) and f 0 (980). The decay width of f 0 (980) is not well determined by experiments, and it is listed as 10−100 MeV. We use the middle values of the Particle Data Group [62], namely Γ σ = 550 MeV between 400 MeV and 700 MeV. As for the decay constants f f0(500) and f f0(980) , no experimental information is available. There are theoretical estimates on f f0(980) by the QCD sum-rule method. However, they assume the qq configuration for f 0 (980) and their decay-constant values seem to be inconsistent with the Belle data on the differential cross section as shown in Sec. IV A. There is no theoretical estimate on f f0(500) as far as we searched, so it is simply terminated or it is considered as one of the parameters in our analysis. These numerical values are summarized in Table I. Next, we define decay constants and distribution amplitudes for the resonances f 0 (500), f 0 (980), and f 2 (1270). In the reaction γ * γ → π 0 π 0 , the matrix elements of a vector current between the vacuum and these meson states are involved in its cross section. First, the matrix element for the tensor meson f 2 (1270) is expressed by the decay constant f f2 and the distribution amplitude Φ f2 (z, µ) as [58,61] where ε (λ) αβ is the polarization vector of f 2 meson [61], and the higher-twist terms are not explicitly written. The distribution amplitude for f 2 is given by the summation of odd Gegenbauer polynomials due to the C-parity as explained in Eq. (78), and it is expressed as [58,61] whereas it is the even polynomials for the pion as shown in Eq. (71).
In the same way, the matrix elements for the scalar mesons f 0 (500) and f 0 (980) are given by where the distribution amplitude is defined by including the decay constant f f0 for a practical purpose, because the combined quantity of f f0 and the amplitude becomes finite even though f f0 itself vanishes. We define the decay constants f f0 andf f0 by the matrix elements for the vector and scalar operators as [63] f Writing the above vector current at the position x as J µ (x) = ψ(x) γ µ ψ(x) = e ip·x J µ (0)e −ip·x and using the equation of motion, we relate the two decay constants as where m q and mq are quark and antiquark masses. In the f 0 -meson case, the masses are equal (mq − m q = 0). Because of the conservation of the vector current or charge-conjugation invariance, the constant f f0 should vanish f f0 = 0. However, the nonlocal matrix element of Eq. (105) does not vanish at finite Q 2 , whereas it vanishes in the scaling limit Q 2 → ∞ as we explain later in Sec. III G. Comparing Eqs. (105) and (106), we obtain the relation For the scalar mesons with m q = mq, the relation (107) can be used to relate the decay constants. Therefore, according to Ref. [63], we may take that the f 0 distribution amplitude is expressed byf f0 and the Gegenbauer polynomials as Then, the normalization of Eq. (108) is satisfied if B 0 is taken as (mq − m q )/m f0 ≡ 1/µ f0 . The integral of the first term is f f0 and those of the subsequent summation terms vanish identically. The first termf f0 /µ f0 = f f0 vanishes for the f 0 meson, so that it is given by where C 3/2 1 (x) = 3x is used.

G. Scale dependence of resonance contributions
There are finite contributions to the γ * γ → π 0 π 0 cross section from f 2 (1270), f 0 (500), and f 0 (980) at small Q 2 . However, as the Q 2 increases, they become smaller and smaller, and they eventually vanish in the scaling limit Q 2 → ∞. The scale dependence of the distribution amplitude is given by the anomalous dimensions γ n and the leading coefficient β 0 = (11C A − 4T R n f )/3 of the β function with C A = N c and T R = 1/2 as [58,61] where C F = (N 2 c − 1)/(2N c ) with the number of colors N c . Here, the meson f indicates f 0 (500), f 0 (980), or f 2 (1270), and the decay constant f f isf f0(500) ,f f0(980) , or f f2(1270) . One could express the scale evolution separately for the decay constant and the distribution amplitude as [61] The leading Gegenbauer polynomial is taken in Eq. (110), and its anomalous dimension is given by γ 1 = 2 C F /3. This finite anomalous dimension indicates that the distribution amplitudes decrease with increasing Q 2 as shown in Eq. (111). From Eq. (112), it is possible to describe the Q 2 evolution separately for the decay constant and the distribution amplitude. However, the overall scale dependence is given by Eq. (111) in any case. The scale dependence is often attributed only to the decay constant [58,61], namely and the distribution amplitude may be normalized in the scale-independent way as so that it becomes as the leading distribution. In Eqs. (104) and (110), the f 2 and f 0 distribution amplitudes are defined with the scale dependence. However, the scale independent expression of Eq. (115) is used in this work, which means to take the B 1 factor as B 1 = 5/3 [61]. This is a consistent description with Eq. (112). However, it may be somewhat confusing, so that one should remember that the distribution amplitude vanishes in the scaling limit Φ f (z, µ) = 0 at µ → ∞, although the scale-independent expression (114) is often used practically.

H. Gravitational form factors for pion
As shown in Eq. (44), the GDAs probe the ++ component of the energy momentum tensor, and it is expressed by the form factors for π 0 as [31,32] Calculating the + components by using the momentum assignments in Eq. (15) and using its relation to the GDAs in Eq. (44), we obtain In this way, if the GDAs are determined from experimental measurements, the gravitational form factors, consequently gravitational radii, are obtained for the pion. Next, we discuss normalizations of the form factors. Using the Legendre polynomial expressed by ζ as P 2 (cos θ) = [−12ζ(1 − ζ) + 3 − β 2 ]/(2β 2 ), we obtain the integral of Eq. (117) as The right-hand-side of this equation should be equal to the sum −4M π 2(q) given in Eq. (19) at W 2 = 4m 2 π , and it leads to the relations in the scaling limit. Quark and antiquark contributions are added to obtain the form factors: Θ n (s = 4m 2 π ) = i=q Θ n,i (s = 4m 2 π ). Then, such sum of the right-hand side of Eq. (122) is q M π 2(q) . Therefore, the normalizations of the form factors become the momentum fraction carried by quarks and antiquarks in the pion: in the scaling limit. The factor q M π 2(q) is written as R π in some articles. Here, the only the quark contributions are discussed, so that the normalization becomes the quark (and antiquark) momentum fraction. However, if the gluon contribution is added, the relation should be Θ 1 (s = 4m 2 π ) = Θ 2 (s = 4m 2 π ) = 1, which indicate A(s = 4m 2 π ) = 1 and B(s = 4m 2 π ) = −1/4 from Eq. (56). Therefore, our timelike form factors are consistent with the works in Refs. [28,31,32]. We should note that these normalizations are satisfied in the scaling limit. However, the Belle measurements are at finite Q 2 with some resonance effects, so that the actual values contain their effects. In fact, as we show later, they are Θ 1 (s = 4m 2 π ) = Θ 2 (s = 4m 2 π ) ∼ 0.7, instead of q M π 2(q) = 0.5, in our GDA analysis.

IV. RESULTS
From these theoretical preparations, we proceed to the actual analysis of experimental data. Here, the Belle data for γ * γ → π 0 π 0 [7] are used for our study. The invariantmass dependent functions are parametrized with the resonance contributions from f 0 (500), f 0 (980) and f 2 (1270), and it is summarized as where the normalization constant N α is given in Eq. (86). The S and D wave terms are expressed by the contributions from the continuum and the resonances as The timelike form factor for the continuum is given by the cut-off parameter Λ and the power of n − 1 in Eq. (94). The factor n is suggested by the constituent counting rule at high energies and it is n = 2 for the pion. Here, f 0 indicates f 0 (500) and f 0 (980). However, the analyzed Belle data are not sensitive to f 0 (980), so that it is not included in our analysis. The up and down quark GDAs are considered in our analysis, and strange and charm quark contributions are neglected. We assigned a parameter α for the z-dependent functional form of the quark GDAs. This z-dependent function enters into the amplitude A ++ in Eq. (65), and then the integral is given in the form of This integral is expressed as the beta function as B(α, α)/(2α + 1), and it plays a role of overall constant to explain the γ * γ → π 0 π 0 data. Next, we explain the S-and D-wave phase shifts used in our analysis. The phase shifts δ 0 and δ 2 in Ref. [24] seem to work below W = 1 GeV. The region of the center-of-mass energy is 0.525 GeV ≤ √ s = W ≤ 2.05 GeV in the Belle data [7]. In order to analyze the Belle data, we use the S-wave and D-wave ππ phase shifts obtained by Bydzovsky, Kaminski, Nazari, and Surovtsev (BKNS) [55]. Their phase shifts are shown in Fig. 11. They proposed a parametrization of the S-and D-wave phase shifts from analysis of the ππ scattering experimental data in the isospin = 0 channel. Since only the difference of S-and D-wave phase shifts matters for explaining the cross section in our analysis of Eqs. (124) and (125), the difference is also shown by the dashed curve. Above the KK threshold at about 1 GeV, the phase difference is, roughly speaking, a slowly varying function of W . We note that the KK channel opens at the threshold energy 2m K + = 0.987354 GeV, so that only the ππ phase shifts may not be sufficient. To be precise, the KK → ππ phase shifts should be introduced together with the kaon GDAs. We may investigate such details step by step. In our GDA analysis, a simple case is considered by introducing a phase ∆δ(W ) for the S wave above the KK threshold, δ 0 (W ) = δ 0 (W ) BKNS +∆δ(W ), δ 2 (W ) = δ 2 (W ) BKNS with expectation that such effects are included in the modified part. In our analysis, we introduce phase parameters in the S-wave as at W > 2m K . The parameters a δ and b δ are determined by the χ 2 analysis.
A. f0(980) contribution A possible complication or ambiguity is how to determine the decay constantsf f0(500) andf f0(980) , whereas the constant f f2 is relatively well evaluated [58,61]. It is because the internal configurations of f 0 (500) and f 0 (980) are not well known. The evaluation off f0 is done for f 0 (980) only by assuming that f 0 is a qq-type meson, namely (uū+dd)/ √ 2 (≡ nn), ss, or mixture of them [63]. On the other hand, it is known that f 0 (980) is likely to be a tetra-quark meson or KK molecule so as to explain the experimental measurements on f 0 (980) → ππ, f 0 (980) → γγ, and φ → f 0 (980)γ [62,64,66]. The theoretical decay constant is not evaluated unfortunately, as far as we are aware, for the tetra-quark or KK configurations for the f 0 (980). Therefore, a realistic numerical estimate would not be possible for f 0 (980) in comparing with experimental data on γ * γ → π 0 π 0 .
Of course, the f 0 (980) may be viewed as a qq state at high energies, whereas it may be a qqqq one at low energies, because they could mix with each other. In fact, there is an indication from the constituent-counting-rule studies on Λ(1405) in comparison with the experimental data on γ + p → Λ(1405)+ K + that Λ(1405) looks pentaquark state (qqqqq) at low energies, whereas it could be an ordinary three-quark one (qqq) at high energies [5]. There is a possibility that the situation could be the same for f 0 (980) on the energy-dependent composition.
In any case, let us simply assume the decay constant f f0(980) by taking the qq-type estimate in the QCD sum rule in order to illustrate the situation and the issue. As we will show later, the optimum value for the parameter α is roughly given by α ∼ 1. We find that it is difficult to accommodate the f 0 (980) resonance with this parameter value. Obtained cross sections are compared with the Belle data at Q 2 = 8.92 GeV 2 and cos θ = 0.1 by taking α = 0.5, 1.0, and 2.0 in Fig. 12. Here, the decay constant f f0(980) (Q 2 = 1 GeV 2 ) = 0.104 GeV was obtained in Ref. [63] by considering the u and d quark contributions to the GDAs withf n = 0.35 GeV at Q 2 = 1 GeV 2 , the mixing angle θ f0 = 32.5 • , which is the middle of 25 • < θ f0 < 40 • , between | nn and | ss , B 1 = −0.92 by the QCD sum rule [63], and the conversion factor 18/30 for the distribution amplitude from Eq. (110) to Eq. (115). The Q 2 evolution is also taken into account by using Eq. (113) from Q 2 = 1 GeV 2 to 8.92 GeV 2 in order to compare with the Belle data at Q 2 = 8.92 GeV 2 . From the comparison with the Belle data in Fig. 12, we find that the f 0 (980) peak structure is not obvious from the data and that they are not consistent with the theoretical predictions as long as α < 2. Although the figure is one of the kinematical point of the Belle measurements, the comparisons with other data also indicate a similar tendency. Here, we should note that the theoretical curves are shown by assuming the qq configuration of f 0 (980) with the QCD sum rule estimate for the decay constant. These results should suggest that f 0 (980) could not be understood mainly by the qq configuration. It is possibly a tetra-quark (or KK molecule) state as widely known in the hadron-physics community. The decay constant could be very small if it is a tetra-quark (qqqq) type because the decay width is proportional to the matrix element of a bilocal operator. Since the data do not show an obvious f 0 (980) peak structure and a theoretical estimate is not available for the decay width by the tetra-quark picture, we do not include f 0 (980) in our numerical analysis in the analysis of Sec. IV B. If the measurements become more accurate in future, one may consider to include this contribution.
The f 0 (980) effect is not conspicuous in the Belle data on the differential cross section, for example, in Fig. 12.
However, it appears in the total cross section [7] and the γ → f 0 (980) transition form factor is investigated by using the Belle data [65]. Such studies indicate that f 0 (980) is consistent with the qq configuration, which is different from our finding in Fig. 12. Because of these conflicting results, the f 0 (980) contribution and its internal configuration are not well understood.

B. Analysis results
The GDAs are expressed by a number of parameters, which are obtained by a χ 2 analysis of Belle experimental measurements on γ * γ → π 0 π 0 . The resonance part is fixed as much as possible by other experimental and theoretical studies, and the used values are listed in Table  I. The large uncertainty comes from the values of the decay constants,f f0(500) ,f f0(980) , and f f2(1270) , especially for the f 0 mesons. There are a number of reliable theoretical studies on f f2(1270) . In Sec. IV A, we explained that the current QCD sum rule estimate forf f0(980) is much different from the Belle measurements if it is assumed as a qq state. There is no available estimate, as far as we are aware, for the decay constant in the tetra-quark picture for f 0 (980). In any case, the data do not show a clear signature of f 0 (980) in the W ∼ 1 GeV region, so that f 0 (980) is not included in the following analysis. Furthermore, there is no theoretical estimate on the decay constant forf f0(500) . We may simply assume that it is same as the f 0 (980) value; however, the results are inconsistent with the Belle data in the same way with the f 0 (980) case. Therefore, we consider two options in our studies: (set 1) Analysis without f 0 (500) and f 0 (980): The GDAs are expressed by the parameters for only the continuum and f 2 (1270), and they are determined by the χ 2 analysis. (set 2) Analysis with f 0 (500) and without f 0 (980): The decay constant f f0(500) is considered as an additional parameter to be determined from the experimental data in addition to the parameters in the set 1. For the decay constants, the Q 2 evolution is taken into account by using Eq. (113) and taking the average scale of the Belle experiment as Q 2 = 16.6 GeV 2 , which is a simple average of the minimum and maximum values, 8.92 and 24.25 GeV 2 , in the analyzed data in this work. By considering the factorization condition of Eq. (12), only the large Q 2 data with Q 2 ≥ 8.92 GeV 2 are used in our analysis. Furthermore, the higher-order and highertwist terms A +− and A 0+ do not contribute significantly at large Q 2 . The Q 2 values to satisfy this condition are Q 2 =8.92, 10.93, 13.37, 17.23, and 24.25 GeV 2 in the Belle measurements. In each Q 2 , the pion angles are cos θ =0.1, 0.3, 0.5, 0.7, and 0.9 as listed in Table II. In each bin of Q 2 and cos θ, there are 22 data points, so that the total number of data is 550.
The GDAs are expressed by the three kinematical variables, z, ζ, and W 2 without the scale Q 2 by considering the scaling region. Actual experiments are done at finite Q 2 , so that the GDAs extracted from the measurements may depend on Q 2 . In order to check the Q 2 dependence of the Belle data, we show the quantity (Q 2 + s)dσ/βd(cos θ) in Figs. 13 and 14 for cos θ = 0.1 and cos θ = 0.5, respectively by choosing W = 0.525, 0.975, and 1.550 GeV. According to Eq. (65), there is no scale dependence for this quantity in the scaling limit. As shown in these figures, the Belle data are not very accurate at this stage to discuss whether Q 2 dependence exists. However, there are tendencies for the scaling within the errors. The Q 2 variations may be seen at Q 2 < 6 GeV 2 at W = 1.55 GeV and cos θ = 0.5 in Fig. 14; however, such data are irrelevant in our analysis because only the data with Q 2 ≥ 8.92 GeV 2 are used. In the analysis 1, there are four parameters, α, Λ, a δ , and b δ , and the others are fixed. For example, n = 2 is taken by the constituent-counting rule, and q M π 2(q) = 0.5 from pion-structure function studies. The f 0 (500) contribution is terminated by takingf f0(500) = 0. The obtained parameter values are listed in Table III. A reasonable fit is obtained in this analysis with χ 2 /d.o.f. = 1.22. Assigning the decay constantf f0(500) as an additional parameter in the analysis 2, we obtained a better agreement with the data with χ 2 /d.o.f. = 1.09. In both cases, the parameter α is close to the asymptotic value α = 1. For the pion distribution amplitude, a more concave functional form is suggested at finite Q 2 [45,53,67]. However, the pion distribution amplitude is related to the C-odd GDAs as shown in Eq. (66), and our current analysis is for the C-even GDAs, so that there is no direct connection.
The cutoff parameter is in the range 1.6 < Λ < 2.0 GeV, which is larger than the cutoff of the nucleon's electromagnetic form factors. The set 2 provides a better description of the Belle data, as indicated by the χ 2 /d.o.f. value, especially at small W (< 0.8 GeV). In both analyses, the values of a δ and b δ stay at almost same values, a δ ≃ 3.8 and b δ ≃ 0.4. In order to explain the Belle data, the decay constant of f 0 (500) isf f0(500) = 0.0183 GeV at Q 2 = 16.6 GeV 2 . It becomesf f0(500) = 0.0246 GeV at Q 2 = 1 GeV 2 , and this value is much smaller than the one forf f0(980) with the qq picture (f f0(980) = 0.104 GeV at Q 2 = 1 GeV 2 ) [63].
The actual comparisons with the Belle data sets are shown in Figs. 15 and 16 for Q 2 =8.92, 13.37, 17.23, and 24.25 GeV 2 and cos θ =0.1 and 0.5. The dashed and solid curves are our theoretical results for set 1 (without f 0 (500)) and set 2 (with f 0 (500)). There is a dip around W = 1 GeV, which is caused by cancellations between the S-and D-wave terms. The f 2 (1270) contribution is obvious at cos θ = 0.1 but it is relatively suppressed at larger cos θ (= 0.5). As mentioned before, the f 0 (980) effects do not appear in the data. However, since the χ 2 value is slightly smaller in the analysis set 2 in comparison with the set-1 value, the f 0 (500) could be needed for interpreting the data in the small W range (W < 0.8 GeV). The whole cross section decreases with increasing Q 2 as shown in Figs.15 and 16. Especially, at reasonably large cos θ, the f 2 resonance effects becomes small. Due to the scale dependence of the decay constantsf f0 and f f2 , the resonance contributions should become small in comparison with the continuum as Q 2 becomes large. At high-energy e + e − colliders such as the international linear collider (ILC), large Q 2 measurements should be done and such experiments are suitable for probing the continuum part of the GDAs. They are valuable for the studies of the GDAs as one of three dimensional structure functions and their relations to the GPDs.
In order to see each term contribution to the γ * γ → π 0 π 0 cross section, we show the cross section solely coming from f 0 (500), GDA continuum, or f 2 (1270) in Fig. 17 by terminating other terms and the phase shifts The dashed and solid curves indicate our analysis results for set 1 (without f0(500)) and set 2 (with f0(500)), respectively. in Eq. (124) for the kinematics of Q 2 = 8.92 GeV 2 and cos θ = 0.1. In the solid curves, the phase shifts are also terminated, whereas the dashed curve indicates the GDA continuum with the phase shifts. For example, the solid GDA continuum curve is obtained by settinḡ f f0(500) = f f2(1270) = 0 and δ 0 = δ 2 = 0. Here, the parameters of the set 2 are used for drawing these curves. In comparison with the solid curve in Fig. 15 for Q 2 = 8.92 GeV 2 and cos θ = 0.1, these distributions seem to be very small. However, the continuum and f 2 contribute to the cross section constructively with almost the same magnitude, so that each contribution is about 1/4 of the cross section of Fig. 15 if other terms are terminated. As expected, f 0 (500) contributes only in the low-energy region of W < 0.8 GeV, and it is much smaller than the continuum according to the set-2 analysis. However, it depends on the f 0 (500) decay constant, which is taken as one of the parameters in our analysis because of the lack of theoretical information. The f 2 (1270) contributes especially in the W = 1.27 GeV region, and its magnitude is comparable to the continuum. The GDA continuum is a slowly varying function of W and it is distributed in the wide W range.

C. Gravitational form factors and radii for pion
Since the optimum GDAs are determined from the Belle data, the timelike gravitational form factors are calculated by Eqs. (119) and (120). Their absolute values are shown in Fig. 18, and individual real and imaginary parts are in Fig. 19. The form factor Θ 2 comes from the D-wave contribution and it is peaked at the f 2 (1270) position, whereas the function Θ 1 has a dip due to the interference between the S-and D-wave terms. The imaginary part of Θ 2 is peaked at the f 2 (1270) resonance and its real part changes the sign. The real and imaginary parts of Θ 1 have both features on the interference and the f 2 (1270) resonance. As for the electric form factor of the pion in the timelike region, there are recent theoretical studies by the holographic QCD and lattice QCD [68]. In order to find the space distributions and gravitational radii, the timelike form factors should be transformed to the spacelike one by using the dispersion relation of Eq. (51). The obtained spacelike form factors are shown in Fig. 20. They are slowly decreasing function of −t and the slope is steeper for Θ 1 than the one for Θ 2 due to the additional S-wave term in Eq. (119).
The substantial difference between the form factors    certainly contradicts to the soft-pion theorem [32,69] which guarantees that Goldstone bosons in gravitational field are insensitive to the scalar curvature [70]. As the gravity is coupled to the conserved energy-momentum tensor including the gluon contributions, it may be the signal that gluon GDA, whose contribution to the considered two-photon process is suppressed, is essential.
Then, the gravitational densities and their radii are calculated by Eqs. (53) and (54), respectively. The gravitational densities ρ 1 (r) and ρ 2 (r), which are obtained from Θ 1 (t) and Θ 2 (t), respectively, are shown for the pion in Fig. 21. It is known that the spacelike electric form factor of the proton is known as the dipole form F p (q) = 1/(1 + q 2 /Λ 2 ) 2 , which leads to the exponential charge density ρ p (r) = (Λ 3 /(8π))e −Λr by the Fourier transform. Typical functional forms of charge densities and form factors are given in Table IV for hadrons and nuclei. The charge form factors and densities of light nuclei are typically given by the Gaussian functional forms, whereas the densities become flat ones for large nuclei. The pion form factor is roughly given by the monopole form F π (q) = 1/(1 + q 2 /Λ 2 ) as suggested by the constituent counting rule, and its space distribution is given by the Yukawa form ρ π (r) = (Λ 2 /(4πr))e −Λr . It is a divergent function as r → 0, so that it is more appropriate to show the density by 4πr 2 ρ(r) rather than ρ(r) itself as usually done for the nucleons and nuclei.
Therefore, the Θ 2 reflects the mass (energy) distribution in the pion. The µν = ij (i, j = 1, 2, 3) components are expressed by the pressure p(r) and shear force s(r) as T ij q ( r ) = p q (r) δ ij + s q (r) Using the definition of the energy-momentum-tensor form factors, we find that p(r) and s(r) are expressed by Θ 1 . Namely, the Θ 1 is the mechanical form factor which contains information on the pressure and shear force. The conservation of the energy-momentum tensor ∂ µ T µν = 0 indicates the stability condition for the pressure p(r) as [29,71] ∞ 0 dr r 2 p(r) = 0.
one B 20 . In physics, the pressure and shear-force distributions have different nature from the mass distribution. We should note that there is uncertainty in our analysis in the sense that only the relative phase δ 0 (W ) − δ 2 (W ) affects the cross section; however, their absolute phases are not. It means that the phase ∆δ(W ) could be attributed to δ 2 instead of δ 0 in Eq. (126). We repeated our χ 2 analysis with this extreme option and obtained the radius values as r 2 mass = 0.56 fm and r 2 mech = 1.56 fm. Therefore, it is fair to state at this stage that the evaluated gravitational radii are in the ranges: It is encouraging that similar radii are obtained in totally different analyses. The mass radius is consistently in the range of the charge radius or slightly smaller, and the mechanical radius is larger. It is interesting to find that there is a possibility that the mass radius is different from the charge one as suggested in Ref. [36]. Lattice QCD calculations on the energy-momentum tensor indicate similar tendency that the mechanical radius is larger than the mass radius [72]. Here, we should note the definition difference from our form factor, namely the factor of −4 in B(t) and Θ 1 (t) as explained below Eq. (123). The actual radii are not shown in the lattice calculations [72]; however, spacelike form factors and radii have similar tendencies with our results. In addition, a theoretical estimate of D-term, which corresponds to Θ 1 in our studies, also shows a similar result [73].
Because the pion GDAs were obtained in this work, it is possible to study their relation to the pion GPDs as explained in Secs. II D and II E. In order to discuss the pion GPDs, we need to find appropriate double distributions from the GDAs and then to calculate the GPDs. Since it is a significant work, we leave it as our future project. In addition, the Belle collaboration has been investigating other hadron-pair production processes including pp. Once such data become available, it is possible to determine nucleonic GDAs in comparison with the GPDs obtained in spacelike reactions.
This kind of studies has a bright prospect in the sense that the Belle collaboration has been analyzing other meson productions γ * γ → hh from the two photon. The experimental errors of Figs. 15 and 16 are dominated by the statistical errors. The KEKB was just upgraded to super-KEKB, so that the errors should be much smaller in the near future. Furthermore, if the ILC is realized, the two-photon cross section γ * γ → hh should be obtained in a very different kinematical region, namely at large Q 2 , and the ILC measurement should be valuable for probing especially the continuum part of the GDAs.
For a long time, the GDAs had been considered as a purely theoretical subject. We showed in this work that it becomes possible to investigate the GDAs experimentally with the appropriate theoretical formalism. This study is merely a starting point. Interesting prospects are waiting for us for investigating gravitational physics for hadrons in the quark-gluon level. For example, the equivalence principle indicates that the anomalous gravitomagnetic moment should vanish in the nucleon [74]. Therefore, the equivalence principle could be tested in the microscopic particle physics by investigating the GPDs and GDAs for the nucleon.

V. SUMMARY
The GDAs are one of three-dimensional structure functions, and they are related to the GPDs by the s-t crossing relation. We analyzed the Belle data of the twophoton cross sections γ * γ → π 0 π 0 for determining the pion GDAs. This work is the first work to obtain the GDAs from the actual experimental data, and our results should be valuable for probing the three-dimensional structure of hadrons, especially for future applications to unstable hadrons including exotic-hadron candidates which cannot be used in fixed-target experiments.
Including the f 0 (500) and f 2 (1270) meson contributions to the cross section, we expressed the pion GDAs by a number of parameters which were determined by analyzing the data. The obtained z-dependence is close to the scaling one (α = 1). If we include f 0 (980) contribution with constants estimated by assuming it as a qq state, theoretical differential cross sections are much larger than the Belle measurements. The f 0 (980) meson was not included in our actual analysis. The GDAs contain the timelike gravitational form factors Θ 1 (s) and Θ 2 (s) of the energy-momentum tensor, and we calculated them from the obtained GDAs. The function Θ 2 (s) is determined only by the D-wave part, whereas both S-and D-waves contribute to Θ 1 (s). Therefore, they have different functional behaviors. This is the first time that the gravitational form factors are obtained from actual experimental measurements.
The timelike gravitational form factors are converted to the spacelike ones by the dispersion relation. Then, the gravitational mass and mechanical densities are shown, and their radii are calculated. We obtained r 2 mass = 0.56 ∼ 0.69 fm and r 2 mech = 1.45 ∼ 1.56 fm from the form factors Θ 2 and Θ 1 , respectively. They indicate that the gravitational mass radius is similar or slightly smaller than the charge radius r 2 charge = 0.672 ± 0.008 fm and that the mechanical radius is larger. Future super-KEKB measurements should improve this situation. We hope that this work will open a new field of gravitational physics in the quark-gluon level.