Determination of the generalized parton distributions through the analysis of the world electron scattering data considering two-photon exchange corrections

We determine the valence generalized parton distributions (GPDs) $ H_v^q $ and $ E_v^q $ with their uncertainties at zero skewness by performing a $ \chi^2 $ analysis of the world electron scattering data considering two-photon exchange corrections. The data include a wide and updated range of the electric and magnetic form factors (FFs) of the proton and neutron. As a result, we find that there are no enough constraints on GPDs $ E_v^q $ from FFs data solely though $ H_v^q $ are well constrained. By including the new data of the charge and magnetic radius of the nucleon in the analysis, we show that they put new constraints on the final GPDs, especially on $ E_v^q $. Moreover, we calculate the gravitational FF $ M_2 $ and the total angular momentum $ J^q $ using the extracted GPDs and compare them with the FFs obtained from the light-cone QCD sum rules (LCSR) and Lattice QCD. We show that our results are interestingly in a good consistency with the pure theoretical predictions.

One advantage of GPDs over PDFs is that they provide quantitative information on both the longitudinal and transverse distributions of partons inside the nucleon. In this way, the structure of the nucleon can be investigated in more details using GPDs which include more degrees of freedom. Actually, GPDs (whether polarized or unpolarized) are the functions of three variables; the fraction of momentum carried by the active quark (x), the square of the momentum transfer in the process (t), and the skewness parameter (ξ), which is a measure of non-forwardness of GPDs.
GPDs contain the extensive information on the hadronic structure. In the forward limit, at zero t and ξ, GPDs are reduced to usual PDFs. The important property of GPDs is that GPDs integrated over x are equal to the corresponding FFs [4]. GPDs are also related to the charge and magnetization distributions. Information on the parton angular momenta can be found from Ji sum rules [3] using H q and E q GPDs. More information on GPDs can be found e.g. in [12,21,22].
The complicated structure of GPDs leads to difficulties on their extraction from a χ 2 analysis of the related experimental data like what was done for the PDFs. As mentioned before, integrals of GPDs over x connect them with corresponding FFs. This relation does not depend on the skewness ξ. This property gives possibility to use the reduced GPDs H q , H q and E q at ξ = 0 to extract them from FFs analyses [23,24]. Following this procedure the analyses of the polarized GPDs H q was done in [25,26].
In present study, we determine unpolarized GPDs H q v and E q v by performing a χ 2 analysis of the world electron scattering data presented in Ref. [27] (YAHL18) where two-photon exchange (TPE) corrections have been incorporated. We also include the data of the charge and magnetic radius of the nucleons in the analysis to investigate their impacts on the extracted GPDs.
The content of the present paper is as follows. In Sec. II, we discuss the phenomenological framework which is used to extract GPDs from data as well as the method for calculating uncertainties. Sec. III is devoted to introduce the experimental data which are included in our analysis. In Sec. IV, we present the results and investigate the goodness of fits. We also compare our results with the corresponding ones from other groups. We summarize our results and conclusions in Sec. V.

II. PHENOMENOLOGICAL FRAMEWORK
As mentioned in the Introduction, GPDs are related to the nucleon FFs. In fact, FFs are first moments of GPDs. For example, at zero skewness, the flavor FFs F q i (i = 1, 2) can be written in terms of the proton valence GPDs H q v and E q v for unpolarized quark of flavor q as On the other hand, the Dirac and Pauli FFs of the nucleon, F 1 and F 2 , are expressed in terms of the flavor FFs F q i as where p and n refer to the proton and neutron, respectively. Therefor, by measuring the nucleon FFs, one can obtain useful information about GPDs. However, the experimental results are typically expressed in terms of the electric and magnetic Sachs form factors, G E and G M , as where j = p, n denotes the type of nucleon. In the above equation, m is the nucleon mass and we have G p E (0) = 1, G n E (0) = 0, and G j M (0) = µ j , where µ j is the magnetic moment of the nucleon. According to DK13 study [24], the valence GPDs H q v (x, t) in Eq. (1), can be expressed in terms of the ordinary valence PDFs q v (x) as where an exponential t behavior is considered. To be more precise, we use ansatz (4) to determine the t-dependences of GPDs with the profile function f q (x) which was introduced in [23] to parameterize the x-dependence of quark distribution in impact parameter space. The profile function f q (x) can have a simple or more flexible form. In this work, following the default analysis of DK13, we use the complex form and take the valence PDFs q v (x) from the ABM11 set [28] at the NLO and scale µ = 2 GeV. The physical motivation behind the profile function (5) has been discussed in details in Ref. [23]. Note also that it leads to a better fit of the data compared with other forms as it has been shown in Refs. [25,26]. The contribution of the strange quark in Eq.
(2) is also neglected as suggested by DK13.
For the valence GPDs E q v (x, t), one can consider the same ansatz of Eq. (4) with a profile function which has similar form to Eq. (5). But, for the forward limit one can not use the usual PDFs in this case. Then, we have with For e q v , we consider a parametrization form similar to DK13 where κ u = 1.67 and κ d = −2.03 are computed from the measured magnetic moments of proton and neutron and the normalization factor N q can be obtained from the fact that An important point about the forward limit of the GPDs and the profile functions is that they can not take arbitrary x dependence due to the fact that the densities for quarks and antiquarks with momentum fraction x at a nominal transverse distance b from the proton center must be different. To be more precise, this requirement implies a positivity condition as follows [24] [e q v (x)] 2 8m 2 ≤ exp(1) where ∆q v (x) are the valence polarized PDFs that we take them from the analysis of the NNPDFpol1.1 [29]. It can be concluded from the above equation that we should have g q (x) < f q (x) to preserve the positivity condition. In fact, as we discuss in Sec. IV, it is very difficult to find a best fit with more flexible profile functions, Eqs. (5) and (7), and distribution e q v (x) in such a way that the positivity is preserved. After describing the phenomenological framework that we use to extract GPDs from data, now it is time to introduce the minimization procedure and the method for calculating uncertainties. In order to determine the unknown parameters of the profile functions (5) and (7) as well as the distribution e q v (x) of Eq. (8), we utilize the standard χ 2 minimization method. To this aim, we use the CERN program library MINUIT [30] and minimize the following function as usual, where summation is performed over all data points included in the analysis. In the above equation, E i is the measured value of the experimental data point i, while T i is the corresponding theoretical estimate. The experimental errors δE i associated with this measurements are calculated from systematic and statistical errors added in quadrature. In order to calculate uncertainties, we use the standard Hessian approach [31] in which the uncertainties of desired quantity F are calculated as where the derivatives are taken with respect to the fitted parameters {η i } with the optimum values {η i }. The Hessian matrix H i,j is calculated by MINUIT and provided at the end of fit procedure. The value of ∆χ 2 determines the confidence region, and we use the standard value ∆χ 2 = 1 in this work.

III. DATA SELECTION
One of the main processes that provide crucial information on GPDs is the elastic electron-nucleon scattering. Actually, by measuring this process one can extract the electric and magnetic FFs of the nucleon (or their ratio) which are related to GPDs as explained in the previous section. One of the methods to extract G E and G M from the unpolarized elastic electron-nucleon scattering is the Rosenbluth separation which provides the separate determination of G E and G M . However, it is believed that the effects of TPE are substantial in this extraction method. There is another method to extract electromagnetic FFs in which the correlation between the polarizations of the beam electron and the proton target (or the scattered proton) is used. An advantage of this method is that it is less sensitive to TPE effects. However, a separate determination of G E and G M is not possible in this method and one can only access to their ratio.
If one considers the single-photon exchange approximation, the cross section of the electron-nucleon scattering can be written in terms of Sachs FFs as [27] where again j = p, n. Here, (dσ/dΩ) Mott is the cross section of the recoil-corrected relativistic point particle (Mott), and , τ are the dimensionless kinematic variables In the above equations, Q 2 = −q 2 is the negative of the momentum transfer squared q 2 to the nucleon, θ is the angle of the final state electron with respect to the incident beam direction, E is the initial electron energy, m is the nucleon mass, and E = E/[1 + (2E/m) sin 2 (θ/2)] is the scattered electron energy. The Born cross section of Eq. (13) should be modified by the radiative corrections as where δ depends on the kinematic variables and includes the vertex, vacuum polarization, and TPE corrections.
In this work, in order to extract GPDs H q v and E q v , we use the data presented in YAHL18 analysis [27]. These data include a wide and updated range of the world electron scattering data off both proton and neutron targets. An important advantage of these data is that the TPE corrections have also been incorporated. We use three datasets of YAHL18, namely, world R p = µ p G p E /G p M polarization, world G n E , and world G n M /µ n G D data where G D = (1 + Q 2 /Λ 2 ) −2 with Λ 2 = 0.71 GeV 2 . These datasets include 69, 38, and 33 data points, respectively. In this way, the total number of data points (N pts. ) that are included in the analysis will be 140. Overall, the data cover the −t range from 0.00973 to 10 GeV 2 .
It should be noted, for the case of unpolarized electron-proton (ep) scattering, the YAHL18 analysis is also included the original cross section data (562 and 658 data points for the world and Mainz cross sections, respectively) which cover the −t range from 0.003841 to 31.2 GeV 2 . An interesting idea is investigation of the impact of original cross section data on the extracted GPDs. This can be done by performing two sets of fits of the YAHL18 data; one considering the cross section and FFs data simultaneously and the other by removing the cross section data and considering just the FFs data. However, as mentioned at the end of Sec. II, it is very difficult to find a best fit that preserves the positivity condition Eq. (10). To be more precise, if one releases the program to find the optimum parameters without considering positivity condition, it leads to unphysical results given that the aim of MINUIT is just minimizing the χ 2 function as far as possible. On the other hand, if one try to implement the positivity condition in the main body of the fit procedure to find the optimum distributions which preserve positivity automatically, the fit does not converge. To overcome this problem one needs to run the program repeatedly as we explain in the next section. Since including the cross section data causes the fit to become time-consuming, we abandon these data and consider only the FFs data.
The other experimental observables related to the electromagnetic FFs that can provide crucial information about the small-t behavior of GPDs are the charge and magnetic radius of the nucleons, r jE and r jM , where again j = p, n stands for the proton and neuron, respectively. Actually, their mean squared are defined by where µ j is the magnetic moment of the nucleon. Note that in these differentiations the terms q v (x)f q (x) and e q v (x)g q (x) are appeared that provide more direct access to the profile functions f q (x) and g q (x). It is also possible to get further information about the e q v (x) distributions. Then, in order to investigate the impact of data of the charge and magnetic radius of the nucleons on the extracted GPDs, we also include them (4 data points) in a new analysis.

IV. RESULTS
In this section, we present the results obtained from the χ 2 analysis of the data introduced in the previous section. To be more precise, we first perform some analyses of the electromagnetic FFs data taken from YAHL18 [27] to construct our base fit. Then, we include also the data of the charge and magnetic radius of the nucleons in the analysis to investigate whether they can put further constraints on the extracted GPDs.
A good approach to find the unknown parameters of the parameterized distributions in any χ 2 analysis of the experimental data is performing a parametrization scan as it is usually used in the global analysis of PDFs [33]. We utilize this procedure to find the optimum values of the parameters and also overcome the positivity problem described at the end of the previous section. The main idea of this procedure is releasing the free parameters step by step and scanning the χ 2 to see how much releasing a parameter affects the value of the χ 2 and also the shape of the distributions obtained. One should continue adding parameters until the change in the value of χ 2 becomes less than unity, ∆χ 2 < 1. So, in this procedure, the parametrization form is obtained systematically. Although such a procedure may seem difficult, since it needs to run the program many times, it leads to optimum results which are physical too.

A. Base fit
As mentioned, we first consider only the YAHL18 data of the electromagnetic FFs, namely R p = µ p G p E /G p M , G n E , and G n M /µ n G D . By analyzing these data one can obtain some base sets of GPDs that can be used (and improved) in the next steps of investigation. In order to find the best fit that preserves the positivity condition, we follow two different scenarios.
• Scenario 1: Since the forward limit of the GPDs E, namely e q v of Eq. (8), plays a crucial role to preserve positivity as pointed out by DK13 [24], we start the parametrization scan by releasing parameters α q and β q of Eq. (8) and parameters α q and A q of the profile function f q (x) of Eq. (5). It should be noted that we set α uv = α dv and fix the parameter α uv by the relation α uv − α dv = 0.1 GeV −2 as suggested by DK13. Moreover, for the scale µ 2 at which the PDFs are chosen in ansatz (4) we also follow DK13 and set it as µ = 2 GeV. As the next step, we fix the parameters β uv and β dv on their optimum values obtained from the previous step and release the other parameters of the profile functions f q (x) and g q (x) step by step. In each step we make sure that the positivity is preserved to a good extent (at least for x < 0.8) as well as the fit is converged. We reject those fits that violate strongly the positivity condition and take the ones with the lowest value of the χ 2 . The procedure is continued until the release of all parameters is checked.
• Scenario 2: Another approach to find the best fit that preserves the positivity condition is that we implement the positivity condition in the main body of the fit program. In this way, the distributions obtained preserve the positivity automatically. However, the problem is that the fit is converged hardly in this case. To solve this problem we confine ourselves to implement condition g q (x) < f q (x) instead of Eq. (10). Then, we use the parametrization scan and follow the same procedure described in Scenario 1. The only difference is that we release parameters β uv and β dv from the beginning to the end of the parametrization scan.
Following Scenario 1, we have found two sets of GPDs that preserve the positivity; one with B dv = 0 and the other with B uv = 0 which are called Set 1 and Set 2, respectively. Note that releasing these parameters (as the last parameter) leads to the positivity violation. Following Scenario 2, we have found another GPD set in which all parameters A q , B q , C q , and D q are free (Set 3). We remind the reader that this set has an another advantage, namely the freedom of the parameters β uv and β dv , that makes the error calculations of the extracted GPDs more real. Another point that should be noted is that for both scenarios (all three sets of GPDs) we considered parameters γ q of Eq. (8) to be equal to zero. In fact, we checked this issue by continuing the parametrization scan and releasing these parameters and found that they can not improve the value of χ 2 significantly. On the other hand, taking them as free parameters leads to the positivity violation. Overall, one can say that the fit is more sensitive to the parameters of the up quark distribution so that releasing them leads to more decrease in the value of χ 2 compared with the corresponding ones of the down quark.
The values of the optimum parameters of the profile functions (5) and (7), and distribution e q v (x) of Eq. (8) are listed in Table I for all three analyses described above. The parameters denoted by an asterisk ( * ) have been fixed as explained. According to the results obtained, both scenarios lead to a large value for parameters β uv and β dv . To be more precise, whether they are fixed after the first step of the parametrization scan on their optimum values or released to the end of the parametrization scan, a large value is obtained for them. This may be strange at the first glance. However, as we will see later, this can be attributed to the fact that there are not enough constrains on e q v (x) from the electromagnetic FFs data solely.  (5) and (7), and distribution e q v (x) of Eq. (8) for the analyses described in Sec. IV A. The parameters denoted by an asterisk ( * ) have been fixed during the fit.

Parameter
Set 1 Set 2 Set 3 Table II contains the results of three analyses of the YAHL18 experimental data for electromagnetic FFs introduced in Sec. III that have been performed following Scenarios 1 and 2 described above. This table includes the list of datasets used in the analysis, along with their related observables. For each dataset, the range of −t which is covered by data and the value of χ 2 divided by the number of data points, χ 2 /N pts. are presented. The last row of the table contains the values of the χ 2 divided by the number of degrees of freedom (χ 2 /d.o.f.). As can be seen, the results obtained for all three analyses Set 1, Set 2, and Set 3 are very similar, especially for Set 1 and Set 2 that have been obtained following same scenario. However, analysis Set 3 has led to the smaller χ 2 . Overall, the values of χ 2 /N pts. for the neutron data are very well that indicates the goodness of fit for these data rather than the proton data of R p . Since in the following we are going to compare our results with the corresponding ones obtained from DK13 analysis, it is appropriate to explain more about the differences and similarities of our analysis with DK13. From a methodological point of view, two analysis are similar to a large extent. Although the parametrization form of the profile functions (5) and (7) as well as the distribution e q v (x) of Eq. (8) are same, we use the parametrization scan, as mentioned before, to find the optimum values of the parameters. In our study, the positivity condition is applied more strictly compared with DK13. Actually, it is either checked step by step during the parametrization scan or applied directly in the main body of the fit program. For the case of data selection, as mentioned in Sec. III, we use the data presented in YAHL18 analysis [27] where the TPE corrections have also been incorporated. Although lots of references are the same, there are also some differences. For the case of proton data, we have used 69 R p data, while DK13 have used 54 data points of both R p and G p M from [34] which can be considered as an older version of YAHL18. For the case of neutron data, we have used 38 and 33 data points of G n E and G n M , respectively, while DK13 have used 36 and 21 data points of G n M and R n , respectively. Overall, the new data points (including updated ones) used in our analysis are 21 data points for R p [35][36][37][38][39], 12 data points for G n M [40][41][42][43], and 3 data points for G n E [44][45][46]. For the case of nucleon's radius data, we have used (see Sec. IV B) 4 data points introduced in Eq. 18, while DK13 have only used r 2 nE . Figure 1 shows a comparison between our results for GPD xH u v (x) obtained from three analyses described above with their uncertainties and the corresponding ones from the analysis of DK13 [24] at four t values, t = 0, −1, −3, −6 GeV 2 . As can be seen from the figure, the results are very similar, especially for Set 1 and Set 2. However, Set 3 predicts smaller distribution. It should be noted that at t = 0 our results are completely consistent with DK13 as expected, since the forward limits of GPDs H q v in all analyses have been taken from ABM11 [28] as mentioned before. It is also worth noting that the maximum position of H u v moves to the larger values of x with −t growing, as expected.
A comparison between our results for GPD xH u v (x) obtained from three analyses described in Sec. IV A and the corresponding one from the analysis of DK13 [24] at four t values shown in panel ( Figure 2 shows the same results as Fig. 1 but for GPD xH d v (x). In this case, Set 1 and Set 2 lead again to very similar results even with −t growing. But, there is significant difference between them and Set 3 especially when −t increases. However, Set 3 and DK13 predict similar distributions at all values of −t. An important point that should be mentioned is that the uncertainties from PDFs in Eq. (4) have not been considered in the error calculations of GPDs H q v in Figs. 1 and 2. Actually, since H q v (x) are related to the valence PDFs q v (x), and on the other hand global analyses of PDFs lead to precise valence parton distributions, it is expected that their uncertainties do not affect significantly the error bands of the valence GPDs H q v (x). The results obtained for GPD xE u v (x) are shown in Fig. 3 and compared again with the corresponding one from DK13 at t = 0, −1, −3, −6 GeV 2 . As can be seen, in this case, Set 1, Set 2, and Set 3 are almost similar, but they differ significantly with DK13. To be more precise, they tend to smaller x rather than DK13. This significant difference can be attributed to the difference in forward limits of E u v , namely e u v (x), as can be concluded from the panel (a) of Fig. 3 that shows the results with t = 0.  Considering all results presented in this section, we conclude that our results for GPDs H q v are in better agreement with DK13 than the corresponding ones for E q v . This difference can be due to both different data included in the analysis and different procedure to find the optimum parameters and final distributions. Note that, as mentioned before, for Set 3 we have implemented the condition g q (x) < f q (x) in the main body of the program so that the set of parameters which lead to the violation of this condition is automatically rejected. Overall, we can say that our results show more respect for the positivity than DK13. In fact, we checked this issue and found that for the case of up quark our results preserve the positivity condition Eq. (10) at all values of x, while for the case of down quark it is violated just for x > 0.8. However, Eq. (10) is violated for both up and down quarks if one considers the results of DK13 for x > 0.9 and x > 0.8, respectively.
Another important conclusion can be drawn from the results obtained in this subsection is that there are not enough constraints on GPDs E q v from the electromagnetic FFs data. In other words, considering such data solely does not lead to the universal GPDs of E q v , since one can obtain different results that all preserve the positivity condition and have the same values of the χ 2 (which means that all of them describe the data well to the same extent).

B. Impact of nucleon's radius data
Having a base set of GPDs in hand, now it is interesting to check how they describe the experimental data of the charge and magnetic radius of the nucleons introduced at the end of Sec. III. To this aim, we have used the final sets of GPDs extracted in the previous subsection as well as DK13 GPDs [24] in the calculation of Eq. (17) and compared the results obtained with the related data in Table. III. As can be seen, DK13 predictions are in better agreement with data rather than our results. In fact, none of our sets of GPDs can predict data within the error, while DK13 predictions of r 2 pM and r 2 nE are inside the error bands of the data. This can be due to the inclusion of r 2 nE data in the DK13 analysis, while our sets are obtained just by the inclusion of the electromagnetic FFs data as described in Sec. IV A. According to the above explanations, in this subsection, we also include the data of the charge and magnetic radius of the nucleon in the analysis to see if they can put more constraints on the extracted GPDs, especially of E q v . To this aim, we follow same scenarios explained at Sec. IV A. The only difference is that we include also parameters β uv and β dv from the beginning to the end of the parametrization scan for the case of scenario 1, since it is expected that more constraints are provided by including the radius data.
Following Scenario 1, we find a set of GPDs with D uv = 0 and γ dv = 0 while all other parameters are free. This set is called Set 4. Note that releasing D uv leads to a strong positivity violation and releasing γ dv does not improve the value of total χ 2 so that it even leads to a little enhancement of the χ 2 /d.o.f. value. Following Scenario 2, we find a set of GPDs with D dv = 0 and γ dv = 0 which is called Set 5. Actually, releasing these parameters do not lead to any improvement in the χ 2 /d.o.f. value. Table IV compares the values of the optimum parameters of the profile functions (5) and (7), and distribution e q v (x) of Eq. (8) for Set 4 and Set 5 described above. The parameters have been fixed during the fit are denoted again by an asterisk ( * ). As can be seen, the values of the parameters β uv and β dv are decreased seriously and become more logical after the inclusion of the radius data in the analysis. The results of these two analyses have also been compared in Table V. It is obvious from this Table that both Set 4 and Set 5 have overall a good description of data and their χ 2 /N pts. are very similar. As before, the goodness of fit for the neutron data is better than the proton data. For the case of charge and magnetic radius data, satisfactory results have been obtained, except for the case of r 2 nM .  (5) and (7), and distribution e q v (x) of Eq. (8) for the analyses described in Sec. IV B. The parameters denoted by an asterisk ( * ) have been fixed during the fit.   Fig. 5 shows a comparison between our results for GPD xH u v (x) obtained from two analyses described in this subsection (Set 4 and Set 5) with their uncertainties and the corresponding ones from the analysis of DK13 at four t values, t = 0, −1, −3, −6 GeV 2 . Here, we have also included the recent results obtained from the Reggeized spectator model (RSM) [47]. Note that, as explained by the authors, their results are valid just for the values of −t less than unity. So, we have not plotted the corresponding ones for t = −3, −6 GeV 2 . As can be seen, our results are very similar to DK13, though they are more suppressed with −t growing compared with DK13. Since the forward limit u v (x) is the same for both cases as clearly seen from panel (a), the difference between our results and DK13 at larger values of −t can be attributed to the difference in the profile functions f uv (x) of Eq. (4). Although the RSM result is compatible with other groups at t = 0, it becomes more different at larger values of −t as it is clear from panel (b). Figure 6 shows the same results as Fig. 5 but for GPD xH d v (x). In this case, our results become a little different rather than DK13 with −t growing that indicates some differences in their profile function f dv (x). Note that DK13 predicts again a larger distribution in magnitude compared with our results. Note also that in this case, the RSM prediction remains in good consistency with others even at larger values of −t.
The results obtained for GPD xE u v (x) are shown in Fig. 7 and compared again with the corresponding ones from DK13 at t = 0, −1, −3, −6 GeV 2 as well as the RSM prediction. As can be seen, at t = 0, Set 4 and DK13 are almost similar, but they differ significantly with Set 5. This indicates that Set 4 and DK13 have similar e u v (x) which tend to the larger values of x rather than e u v (x) of Set 5. However, they become also different by increasing the absolute value of t that indicates some differences in their profile function g uv (x) of Eq. (7). The RSM result is more compatible with Set 4 and DK13, though its peak tends to the larger values of x with −t growing. Figure 8 shows the same results as Fig. 7 but for GPD xE d v (x). In this case, our results are in good agreement with each other and some differences appear just with −t growing. However, they differ significantly from DK13 at all values of −t. Actually, our results tend to larger x values and have also larger magnitude compared with DK13. At t = 0, the RSM behaves as our results at small x values, but it becomes different at medium and large values of x. Although, the RSM and DK13 are not in good consistency at t = 0, they become similar at t = −1 GeV 2 .
Comparing the results obtained in this subsection with the corresponding ones from Sec. IV A, one can conclude that the inclusion of the charge and magnetic radius data in the analysis put new constraints one GPDs, especially of E q v . Actually, it is obvious from the results obtained that by inclusion of these data in the analysis, xE q v distributions have significantly shifted to the large values of x, while without these data they tend to localize at medium x. Moreover, xH q v distributions are decreased in magnitude by inclusion the radius data. Note that, comparing with GPDs H q v , there are still some differences between different sets of GPDs of E q v . This indicates that it is necessary to include more precise experimental data in the analysis, especially those can put more constraints on the GPDs E q v .

C. Comparison with other quantities
After extracting different sets of GPDs utilizing different scenarios and including different types of experimental data, it is interesting now to investigate how they describe the other physical quantities which are related to GPDs at zero skewness (ξ = 0).
It is well established now that nth Mellin moment of the GPDs H and E can be expressed as polynomials in ξ of order n + 1 considering the polynomiality property of the GPDs [48]. For example, for the second Mellin moment of the GPD H we have [49] where M 2 (t) and d 1 (t) are the gravitational FFs of the energy-momentum tensor. It is known that M 2 (t) can provide information on the momentum fractions carried by the constituent quarks of the nucleon, while the information on the distribution of pressure and tensor forces inside hadrons can be accessed from d 1 (t) (called the D-term). For more information on gravitational FFs, one may refer to Refs. [50,51] as examples. From Eq. (19), it is obvious that M 2 (t) is related to the GPDs H q (x, t) at zero skewness. Then, it may be an interesting idea to calculate M 2 (t) using our GPDs obtained in Sec. IV B as a function of −t. Note that if both t and ξ are set to zero, the M 2 will be a simple sum of the momentum fractions carried by quarks, since GPDs H q (x) turn into the ordinary PDFs q(x). Although there are no experimental measurements of M 2 (t) unlike the D-term d 1 (t), it is interesting to compare our results with the corresponding ones obtained from the light-cone QCD sum rules (LCSR) [52]. Such a comparison has been shown in Fig. 9 where we have calculated the left hand of Eq. (19) using our GPDs obtained in Sec. IV B (Set 4 and Set 5) and compared them with the results of M 2 (t) obtained using LCSR. It should be noted that since we have extracted only the valence GPDs from the experimental data and our results do not include the contributions from the sea quarks, we have considered an assumption to calculate M 2 (t) that is using the profile functions of the valence quarks for the sea quarks too. Moreover, the LCSR result shown in Fig. 9 is the averaged results of LCSR presented in Ref. [52] that has been evolved to µ = 2 GeV using the the renormalization group equations extended to include mass renormalization (the original results are belonging to µ = 1 GeV).
As can be seen from Fig. 9, the results are in good agreement with each other especially if one considers the uncertainties. Actually, it is very interesting that the result of a phenomenological approach in which GPDs are constrained from the experimental data is in a good consistency with the pure theoretical calculations. This indicates, on the other hand, the validity of the LCSR framework to study the hadrons structure and their properties. Note that the excellent agreement between Set 4 and Set 5 is a result of good constraints on GPDs H q . Figure 10 shows same comparison as Fig. 9 but in the interval 0 < −t < 1.3 GeV 2 and including also the Lattice results taken from Table 15 of Ref. [53]. As can be seen, the consistency is also good in this case considering this fact that there is no explicit results from LCSR for −t < 1 GeV 2 and the extrapolation has been used to extend the calculation to this interval. Another quantity that is related to GPDs is the total angular momentum carried by the quarks inside the nucleon. In the most general case, it can be expressed in terms of GPDs H q and E q at zero skewness for each desired quark flavor q as It should be noted that this relation turns into the famous Ji's sum rule [4,54] at t = 0. Comparing with M 2 , J q may be more interesting for us since it contains also GPD E q . Since our analyses just includes the valence GPDs, we can calculate J u v and J d v . Figure 11 shows a comparison between our predictions for J u v and J d v at limit t = 0 and µ = 2 GeV obtained again using Set 4 and Set 5 and the corresponding ones from PRC88 [55], LHPC [56], Thomas [57], TMD [58], and DK13 [24]. According to this figure, our results for J u v are in more agreement with LHPC, TMD, and DK13 while they have significant difference with PRC88 and Thomas. For the case of J d v , our results are in good agreement with PRC88 and TMD while differ significantly with other groups. Overall, one can conclude from Fig. 11 that our results, especially Set 4, are more consistent with LHPC, TMD and DK13.
Other than Ji's sum rule for each specific flavor q at t = 0, we can calculate Eq. (20) as a function of t and also sum over all flavors to get the total angular momentum of the proton J q . However, since our analyses dose not include contributions from sea quarks we can just calculate J q using the valence GPDs. Figure 12 shows the results obtained using Set 4 and Set 5 and compares them with the corresponding ones from LCSR similar to Fig. 9 for M 2 (t). Overall there is a good agreement between our results and LCSR again. However, the results are more different at small values of −t compared with Fig. 9. Note that the difference between Set 4 and Set 5 in this case comes from the difference between GPDs E q of these two sets. However, Set 4 is more consistent with LCSR. Moreover, some of the differences between our results and LCSR can be attributed to the absence of the sea quark contributions in our calculations.
We have compared the results obtained for J q of the proton at the interval 0 < −t < 1.3 GeV 2 in Fig. 13. As before, the figure includes also the lattice results taken from Table 15 of Ref. [53]. Although the agreement between the results are acceptable considering uncertainties, but the differences are more considerable compared with Fig. 10 for M 2 .

V. SUMMARY AND CONCLUSION
The 3D hadron structure can be accessed through GPDs [1][2][3][4][5] which are measurable in hard exclusive scattering processes. In this work, following the recent works [25,26] performed to determine the polarized GPDs H q v , we determined the unpolarized valence GPDs H q v and E q v with their uncertainties at zero skewness (ξ = 0) by performing some χ 2 analyses of the related experimental data. To this end, we first considered the world electron scattering data presented in Ref. [27] (YAHL18) where TPE corrections have also been incorporated. These data include the world R p = µ p G p E /G p M polarization, world G n E , and world G n M /µ n G D data with overall 140 data points. They cover the −t range from 0.00973 to 10 GeV 2 . Utilizing two different scenarios where the optimum values of the parameters are obtained through a parametrization scan procedure [33], we extracted three different sets of GPDs that all preserve the positivity condition and lead to an acceptable quality of the fit, especially for the neutron data. We compared our GPDs with the corresponding ones obtained from DK13 [24] analysis. Although our results for H q v are in good consistencies with DK13, the results for E q v are different. Overall, we found that there are not enough constraints on FIG. 9. A comparison between our results obtained for M2(t) using Set 4 and Set 5 of GPDs and the corresponding one calculated using LCSR [52].
GPDs E q v from FFs data solely though H q v are well constrained. As the next step, we included the data of the charge and magnetic radius of the nucleons in the analysis to investigate their impacts on the extracted GPDs. Utilizing two different scenarios again, we obtained two final sets of GPDs, namely Set 4 and Set 5, which are in more consistent with DK13, especially Set 4. We shown that the radius data put new constraints one the final GPDs, especially of E q v .
To be more precise, by inclusion of these data in the analysis, xH q v distributions are decreased in magnitude and xE q v are significantly shifted to the large values of x rather than before.
As the final step, we calculated the gravitational FF M 2 (t) and the total angular momentum of proton J q (t) using our final sets of GPDs Set 4 and Set 5 and compared them with the results obtained from LCSR [52] and Lattice QCD [53]. We shown that our results are interestingly in a good consistency with the pure theoretical predictions, especially Set 4. Overall, the differences are more considerable for the case of J q (t) that includes also GPDs E q  Fig. 12, but in the interval 0 < −t < 1.3 GeV 2 and including also the lattice results taken from Table 15 of Ref. [53].
compared to M 2 (t). In addition, we calculated the total angular momentum carried by the valence quarks at t = 0, namely Ji's sum rule [4], and compared our results with the corresponding ones obtained from other groups. Overall, our results, especially Set 4, are more consistent with LHPC [56], TMD [58] and DK13 [24]. According to the results obtained in this paper, we emphasize that in order to get more universal GPDs, especially of E q v , it is necessary to include more precise experimental data in the analysis. In this regard, the future programs like one that will be done at Jefferson Lab [59] can shed new lights on the determination of GPDs.