Constraints on charm-anticharm asymmetry in the nucleon from lattice QCD

We present the first lattice QCD calculation of the charm quark contribution to the nucleon electromagnetic form factors $G^c_{E,M}(Q^2)$ in the momentum transfer range $0\leq Q^2 \leq 1.4$ $\rm GeV^2$. The quark mass dependence, finite lattice spacing and volume corrections are taken into account simultaneously based on the calculation on three gauge ensembles including one at the physical pion mass. The nonzero value of the charm magnetic moment $\mu^c_M=-0.00127(38)_{\rm stat}(5)_{\rm sys}$, as well as the Pauli form factor, reflects a nontrivial role of the charm sea in the nucleon spin structure. The nonzero $G^c_{E}(Q^2)$ indicates the existence of a nonvanishing asymmetric charm-anticharm sea in the nucleon. Performing a nonperturbative analysis based on holographic QCD and the generalized Veneziano model, we study the constraints on the $[c(x)-\bar{c}(x)]$ distribution from the lattice QCD results presented here. Our results provide complementary information and motivation for more detailed studies of physical observables that are sensitive to intrinsic charm and for future global analyses of parton distributions including asymmetric charm-anticharm distribution.


Introduction
The charm-anticharm sea in the nucleon has received great interest in nuclear and particle physics for its particular significance in understanding high energy reactions associated with charm production. Quantum Chromodynamics (QCD), the underlying theory of the strong interaction with quarks and gluons as the fundamental degrees of freedom, allows heavy quarks in the nucleon-sea to have both perturbative "extrinsic" and nonperturbative "intrinsic" origins. The extrinsic sea arises from gluon splitting triggered by a probe in the reaction. It can be calculated order-by-order in perturbation theory if the probe is hard. The intrinsic sea is encoded in the nucleon wave functions.
The existence of nonperturbative intrinsic charm (IC) was originally proposed in the BHPS model [1] and in the subsequent calculations [2,3,4] following the original proposal [1]. Proper knowledge of the existence of IC and an estimate of its magnitude will elucidate some fundamental aspects of nonperturbative QCD. Therefore, the main goal of this article is to investigate the existence of a nonzero "intrinsic" charm of nonperturbative origin in the nucleon. In the case of light-front (LF) Hamiltonian theory, the intrinsic heavy quarks of the proton are associated with higher Fock states such as |uudQQ in the hadronic eigenstate of the LF Hamiltonian; this implies that the heavy quarks are multiconnected to the valence quarks. The probability for the heavy-quark Fock states scales as 1/m 2 Q in non-Abelian QCD. Since the LF wavefunction is maximal at minimum off-shell invariant mass; i.e., at equal rapidity, the intrinsic heavy quarks carry large momentum fraction x Q . A key characteristic is different momentum and spin distributions for the intrinsic Q andQ in the nucleon; for example the charm-anticharm asymmetry, since the comoving quarks are sensitive to the global quantum numbers of the nucleon [5].
IC was also proposed in meson-baryon fluctuation models [6,7]. Because of possible direct and indirect relevance of IC in several physical processes (discussed below), there have been many phenomenological calculations involving the existence of a non-zero IC component to explain anomalies in the experimental data and possible signatures of IC in upcoming experiments [5]. Unfortunately, the normalization of the |uudcc intrinsic charm Fock component in the lightfront wavefunctions (LFWF) is unknown. Also, the probability to find a two-body stateD 0 (uc)Λ + c (udc) in the proton within the meson-baryon fluctuation models cannot be determined without additional assumptions: precise constraints from future experiments and/or firstprinciples calculations are required.
The effect of whether the IC parton distribution is either included or excluded in the determinations of charm parton distribution functions (PDFs) can induce changes in other parton distributions through the momentum sum rule, which can indirectly affect the analyses of various physical processes that depend on the input of various PDFs. On the experimental side, an estimate of intrinsic charm (c) and anticharm (c) distributions can provide important information to the understanding of charm quark production in deep inelastic lp → l cX scattering in the EMC experiment [8]. The enhancement of charm distribution in the measurement the charm quark structure function F c 2 compared to the expectation from the gluon splitting mechanism in the EMC experimental data has been interpreted as evidence for nonzero IC in several calculations [2,3,9,10]. A precise determination of charm and anticharm PDFs by considering both the perturbative and nonperturbative contributions is important in understanding charmonia and open charm productions, such as the J/ψ production at large momentum from pA collisions at CERN [11], from πA collisions at FNAL [12], from pp collisions at LHC [13], and charmed hadron or jet production from pp collisions at ISR, FNAL, and LHC [13,14,15,16]. LHC measurements associated with cross section of inclusive production of Higgs, Z, W bosons via gluon-gluon fusion, and productions of charm jet and Z 0 [17,18,19,20], J/ψ and D 0 mesons at LHCb experiment [13] can also be sensitive to the intrinsic charm distribution. The J/ψ photo-or electroproductions near the charm threshold is believed to be sensitive to the trace anomaly component of the proton mass, and some experiments have been proposed at JLab [21] as well as for the future EIC to measure the production cross section near the threshold. The existence of IC in the proton will provide additional production channels and thus enhance the cross section, especially near the threshold. Similarly, open charm pro-duction will also be enhanced by IC. If c andc quarks have different distributions in the proton, the enhancements on D andD productions will appear at slightly different kinematics. IC has also been proposed to have an impact on estimating the astrophysical neutrino flux observed at the IceCube experiment [22].
In global analyses of PDFs, one parametrizes the distributions at an initial scale µ 0 and evolves them to higher scales to fit data. There are different approaches to deal with heavy quarks in which a transition of the number of active quark flavors is made at some scale around the charm quark mass µ c ∼ m c [23,24,25,26]. The transition scale defines where the extrinsic charmanticharm sea enters. However, the intrinsic charmanticharm sea can exist even at a lower scale. In many global fits [27,28,29,30,31,32] the charm quark PDF is set to zero at µ c , but this is an assumption of no IC. In recent years, several PDF analyses started to investigate the possibility of nonzero charm and anticharm distributions at the scale µ c [33,34,35,36,37], but none of them can provide conclusive evidence or exclusion for the intrinsic charm due to the absence of precise data. We need to point out that a nonzero charm quark PDF at µ c is not necessarily evidence of intrinsic charm, because such estimation depends on the heavy-quark scheme and the µ c value used in the fit. This explains the speculation in [33] that the estimation of intrinsic charm may strongly depend on the choice of the transition scale µ c . Fortunately, there is an ideal quantity, the asymmetric charm-anticharm distribution [c(x) −c(x)], which would be a clear signal for IC. Such asymmetry is allowed in QCD, since the interaction of quark-quark is different from that of quark-antiquark. The nucleon has nonzero quark number, the number of quarks minus the number of antiquarks, and thus the c quark and thē c quark in a nucleon would "feel" different interactions, leading to an asymmetric charm-anticharm distribution. Although the absence of such asymmetry does not exclude intrinsic charm, a nonzero [c(x) −c(x)] can serve as a strong evidence, because the extrinsic part of such asymmetry arising at the next-to-next-to-leading order level is negligible [38].
Although the global fits in [33,34,35,36,37] consider the possibility of IC, all these fits assume [c(x) − c(x)] = 0; constraints on [c(x) −c(x)] have been warranted in the NNPDF global fit [36]. It was found in the global fit [36] that a precise and accurate parametrization of the charm PDFs will be useful for more reliable phenomenology using the data from LHC experiments and will eliminate possible sources of bias arising from the assumptions of only perturbatively-generated charm PDFs. It is therefore important to determine the magni-tude of the IC, small or large, namely if [c(x)−c(x)] 0, and how or whether the IC will have significant effect in the physical processes [11,12,13,14,15,16,17,18,19,20]. A precise and accurate knowledge of the IC will also have a direct impact on determining the unknown normalization constants of different model calculations associated with the nonperturbative c(x) andc(x) distributions.
An important question to ask is whether the firstprinciples lattice QCD (LQCD) calculation can provide some constraints or complementary information regarding the existence of IC. Recently, LQCD has provided several calculations of hadron structure complementary to experiments. For example, there have been recent LQCD calculations of the strange (s) quark electromagnetic form factors [39,40,41,42,43,44] with higher precision and accuracy than previously attained by experiments. These LQCD calculations provided indirect evidence for a nonzero strange quarkantiquark asymmetry in the nucleon by pinning down the nonzero value of the strange electric form factor G s E (Q 2 ) at Q 2 > 0. The determination of G s E,M (Q 2 ) from the lattice QCD calculation in [40] has led to precise determination of neutral current weak axial and electromagnetic form factors [42,45]. Using LQCD results from [40,41] as constraints, it has been shown recently in [46], within the light-front holographic QCD (LFHQCD) approach [47,48,49,50] and the generalized Veneziano model [51,52,53], that the [s(x) −s(x)] distribution can be computed within this nonperturbative framework. The results from [46] lead to a negative value of [s(x) −s(x)] at small-x and positive at largex, showing the possibility of applying LQCD results for the phenomenological study of the [s(x) −s(x)] asymmetry in the absence of precise experimental data and global fits of the strange quark PDFs.
The main goal of this article is to calculate the charm electric and magnetic form factors, G c E (Q 2 ) and G c M (Q 2 ), from LQCD at nonzero momentum transfer and discuss their connection to the existence of IC and a nonzero [c(x) −c(x)] asymmetry distribution in the nucleon. Using the LQCD calculation of G c E,M (Q 2 ) as constraint, we determine the [c(x) −c(x)] distribution using the nonperturbative framework described in [54]. We note that the electromagnetic current is odd under charge conjugation and the Dirac form factor F c 1 (Q 2 > 0) provides a measure of the c-quark minus thec-quark contribution due to the opposite charges of the quark and antiquark. While F c 1 (Q 2 = 0) = 0, according to the requirement that an equal number of charm and anticharm quarks are required by the quantum numbers of the nucleon, a positive F c 1 (Q 2 ) at Q 2 > 0 implies that the c-quark distribution is more centralized than thec quark distribution in coordinate space. This, in turn, results in a [c(x) −c(x)] asymmetry in momentum space, thereby providing possible evidence for the nonperturbative IC in the nucleon. On the other hand, a nonzero charm Pauli form factor F c 2 (Q 2 ) and a nonzero charm magnetic moment µ c M 0 are consequences of a nonzero orbital angular momentum contribution to the nucleon from charm quarks [55]. This can be understood in the inherently relativistic LF formalism where a nonzero anomalous magnetic moment requires to have orbital angular momentum L z = 0 and L z = 1 Fockstates components in the LFWF.
The remainder of this article is organized as follows: In Section 2, we briefly discuss the LQCD calculation of the charm electromagnetic form factors G c E,M (Q 2 ) and determine these in the physical limit. In Section 3, we use LQCD results for G c E,M (Q 2 ) as input for quantitative analysis of the [c(x) −c(x)] asymmetry distribution in the nucleon within the specific framework of LFHQCD and the generalized Veneziano model. We also present a brief qualitative discussion of our results in Section 4.

Lattice QCD calculation of G c E,M (Q 2 )
We present in this section the first lattice QCD calculation of the charm quark electromagnetic form factors in the nucleon. This first-principles analysis requires a disconnected insertion calculation. By "disconnected insertion," one refers to the nucleon matrix elements involving self-contracted quark graphs (loops), which are correlated with the valence quarks in the nucleon propagator by the fluctuating background gauge fields (n.b., notice the difference with the usual term "disconnected diagram" used in the continuum Quantum Field Theory literature). Numerical expense and complexity of the disconnected insertion calculations in LQCD, the deficit of good signal-to-noise ratio in the matrix elements, and the possibility for a very small magnitude of the c quark matrix elements make it difficult to obtain a precise determination of G c E,M (Q 2 ). We, therefore, need to accept several limitations while performing this calculation. For example, the data is almost twice as noisy compared to the matrix elements of the strange electromagnetic form factors G s E,M (Q 2 ) [40, 41, 42] and we do not see any signal for one of the gauge ensembles (32ID with a lattice spacing of a = 0.143 fm [56]) used in the previous calculations [40,42]. Moreover, we are only able to perform the widely used two-states summed ratio fit of the nucleon three-point (3pt) to two-point (2pt) correlation functions instead of a simultaneous fit to the summed ratio and conventional 3pt/2pt-ratio as was done in [42]. The reason is that the 3pt/2pt-ratio fit for extracting G c E,M (Q 2 ) is not stable for the ensemble at the physical pion mass m π = 139 MeV. We also keep in mind that the O(m 2 c a 2 ) errors associated with the lattice spacing can be larger than the case for the G s E,M (Q 2 ) matrix elements.
Our calculation comprises numerical computation with valence overlap fermion on three RBC/UKQCD domain-wall fermion gauge configurations. We use 17 valence quark masses in total, corresponding to pion masses in the range m π ∈[139,330] MeV across these ensembles to explore the quark-mass dependence of the charm electromagnetic form factors. Parameters for each gauge ensemble used in this work are given in Table 1. The details of the numerical setup of this calculation can be found in [58,59,60,42]. The charm quark mass m c was determined in a global fit on the lattice ensembles with β = 2.13 fm and 2.25 using inputs from three physical quantities, such as M D * s , M D * s − M D s , and M J/ψ in [61]. Our statistics are from approximately 100k to 500k measurements across the 24I to 48I ensembles. Other than elaborating on the numerical techniques, we refer the readers to previous work [42] for a detailed discussion of the similar numerical techniques which have been used for this calculation.
We extract G c E,M (Q 2 ) matrix elements using the twostates fit of the nucleon 3pt/2pt summed ratio SR(t 2 ): Here, R(t 2 , t 1 ) is the 3pt/2pt-ratio, t 0 and t 2 are the source and sink temporal positions, respectively, and t 1 is the time at which the bilinear operator V µ (x) = c(x)γ µ c(x) is inserted, t and t are the number of time slices we drop at the source and sink sides, respectively, and we choose t = t = 1. C i are the spectral weights involving the excited-state contamination. Ideally, ∆m is the energy difference between the first excited state and the ground state but in practice, this is an average of the mass difference between the proton and the lowest few excited states. The excited states in the SR(t 2 ) fit fall off faster as e −∆mt 2 compared to two-states fit case of the 3pt/2pt-ratio where the excited-state falls off at a slower rate as e −∆m(t 2 −t 1 ) . These faster-decreasing excited-state effects allow for fitting the matrix elements starting from short time extents, as was demonstrated in [62]. We present the matrix elements of G c E,M (Q 2 ) obtained from the fit Eq. (1) in Fig. 1 and Fig. 2.   Ensemble With the extracted 96 matrix elements from three gauge ensembles (separately for G c E (Q 2 ) and G c M (Q 2 )) at different pion masses and Q 2 , we perform a simultaneous model-independent z-expansion fit [63,64] to the G c E,M (Q 2 ) in the momentum transfer range of 0 ≤ Q 2 ≤ 1.4 GeV 2 and perform chiral, continuum (lattice spacing a → 0), and infinite volume limit (lattice spatial extent L → ∞) extrapolation to obtain the form factors in the physical limit. For such a fit to the G c E (Q 2 ), we adopt the following fit form where In fit Eq. (2), m π,vs is the partially quenched pion mass m 2 π,vs = 1/2(m 2 π + m 2 π,ss ) with m π,ss the pion mass corresponding to the sea quark mass, and m J/ψ masses for the lattice ensembles are obtained in [61] and extrapolated to the physical value m J/ψ = 3.097 GeV [65]. A 4 includes the mixed-action parameter ∆ mix [66]. The volume correction in fit (2) has been adopted from [67] to best describe the data. We use t cut = m 2 J/ψ , the pole of cc pair production. We note that this choice is different from the fit to the strange quark form factor where the t cut is chosen at 4m 2 K , because the mass of two kaons is less than the mass of φ, while the mass of two D mesons is greater than the mass of J/ψ. One may also consider η c , which is a bit lighter, but J/ψ is more likely to be produced from a vector current.
Inclusion of higher-order terms beyond k max = 4 has no statistical significance and are not considered in the z-expansion fit (2). We obtain χ 2 /d.o.f. = 1.17 for the fit (2) and the fit parameters are The significant increase of the uncertainty in the physical value of G c E (Q 2 ) at larger values of Q 2 is due to the fact that the data points on the 24I and 32I ensembles are at much heavier pion mass compared to the matrix element at the physical m π = 139 MeV on the 48I ensemble and there exist no LQCD data points at Q 2 ≥ 0.31 GeV 2 on the 48I ensemble. We also see a similar feature for the G c M (Q 2 ) shown in Fig. 2. The cyan band in Fig. 1 represents G c E (Q 2 )| physical in the physical limit after the quark mass, finite lattice spacing and volume corrections have been implemented using the fit parameters listed above.
The systematic uncertainty is estimated by calculating the differences between the G c E (Q 2 )| physical and G c E (Q 2 ) obtained from the fit Eq. (2) by considering the corrections of the A 3 , A 4 , A 5 terms from the m J/ψ -value on the 24I ensemble obtained in [61], the smallest lattice spacing from 32I ensemble, and the 48I ensemble with the largest volume, respectively.
To obtain the charm quark magnetic form factor G c M (Q 2 ) in the physical limit and in the 0 ≤ Q 2 ≤ 1.4 GeV 2 momentum transfer region, we adopt the following empirical fit form with the volume correction term adopted from [68]: Similar to the z-expansion fit to the charm electric form factor, we limit the k max in our fit to the value for which an additional term in the z-expansion has no statistical significance and obtain k max = 3. With the χ 2 /d.o.f. = 1.14 in the fit (4), we obtain the following fit parameters λ 0 λ 1 λ 2 −0.00127 (38) 0.054 (20) −0.86(70) The systematic uncertainty of G c M (Q 2 )| physical is obtained in a similar way as for the case of G c E (Q 2 )| physical . The electric and magnetic radii can be extracted from the slope of the G c E,M (Q 2 ) form factors as Q 2 → 0: Using the λ 1 -values from above, we obtain While the existence of a nonzero G c E (Q 2 ) (and thus the Dirac F c 1 (Q 2 ) form factor) is related to a nonzero asymmetry of the [c(x) −c(x)] distribution, a nonzero G c M (Q 2 ) (and thus Pauli F c 2 (Q 2 ) form factor) immediately implies that there is a nonzero orbital angular momentum contribution to the nucleon from the charm quarks. The Pauli form factor F c 2 (Q 2 ) has the LFWF representation as the overlap of the states differing by one unit of orbital angular momentum. It is thus closely related to the spin sector of the charm quark sea in the proton. F c 2 (Q 2 ) is also given by the first x moment of the generalized parton distributions E c (x, ξ, t), t = −Q 2 , which contributes to the second term of Ji's sum rule [69]. At Q 2 = 0, F c 2 (0) is equal to the magnetic moment µ c M since F c 1 (0) = 0 is fixed by the quark number. Therefore, our result µ c M = 0.00127(38) stat (5) sys indicates a nontrivial role of the charm quark sea in understanding the spin content of the proton.

A nonperturbative model for computing the intrinsic [c(x) −c(x)] asymmetry in the nucleon
The existence of IC in parton distributions was first proposed in the BHPS model [1] and then followed by various phenomenological models. Although a modelindependent determination of IC distributions from the charm quark form factors is not possible at the moment, the underlying physics, governed by QCD, does imply some connections among them, such as the QCD inclusive-exclusive connection [70,71,72]: It relates hadron form factors at large Q 2 to the hadron structure function at x → 1, leading to simple counting rules. Furthermore, since the charm and anticharm carry opposite charges, the charm quark form factor measures the difference of the transverse charge density [73] between the charm quark and the anticharm quark. For a positive charge form factor, like those shown in Fig. 1, the charm quark distribution is more spread than the anticharm distribution in q T -space, where q T is the Fourier conjugate variable of the transverse coordinate b T and q 2 T = Q 2 . As a feature of Fourier transform, the charm quark density ρ(b T ) is more centralized in b T -space than the anticharm quark. Representing the density as the square of the wave function, ρ(b T ) = |ψ(b T )| 2 , one can easily find the k T -space wave function ψ(k T ) is more spread out for the charm quark, where k T is the intrinsic transverse momentum. Thus the charm quark favors carrying higher momentum than the anticharm quark for a positive charge form factor. As a result, the charmanticharm asymmetric distribution function [c(x) −c(x)] will favor negative values in the low-x region and positive values in the high-x region. We note that a strict definition of the transverse charge density is given by the Fourier transform of the Dirac form factor F c 1 (Q 2 ), which dominates the charge form factor G c E (Q 2 ) in the low-Q 2 regime.
For a quantitative estimation of the charm-anticharm asymmetric distribution function [c(x) −c(x)], one currently has to rely on additional assumptions, although the form factor result does indicate some qualitative features of the distribution function based on the discussions above. Here we take the nonperturbative phenomenological model in [46], which relates the form factor and parton distribution functions with minimal parameters. The formalism is based on the gauge/gravity correspondence [74], light-front holographic mapping [49,54,50], and the generalized Veneziano model [51,52,53]. In the following, we refer to this model shortly as LFHQCD. The charm quark Dirac and Pauli form factors are given by [75] with where τ is the number of constituents of the Fock state component, α(t) is the Regge trajectory, and N τ is a normalization factor. The light front holographic approach leading to these results is based on the underlying conformal algebra which leads to linear trajectories. The spin-flavor coefficients c τ and χ τ are parameters to be determined from the LQCD computation of G c E (Q 2 )| physical and G c M (Q 2 )| physical to obtain F c 1 (Q 2 ). The constraint that the numbers of charm and anticharm quarks are identical for each Fock state component has been implemented in Eq. (7).
Using the integral representation of the Beta function B(τ − 1, 1 − α(t)) in Eq. (9), the form factor is expressed in a reparametrization invariant form [54] where w(x) is a flavor independent function with w(0) = 0, w(1) = 1 and w (x) ≥ 0. Then the asymmetric charmanticharm distribution function is where We should note here that the coefficients c τ in Eq. (11) are the same as those in Eq. (7). Therefore, once they are determined by the form factor, one can make predictions for the distribution functions. The form factors and distribution functions above are derived at the massless quark limit and one may have different approaches to incorporate quark mass corrections. For small quark masses (up, down and strange) the latter can be treated perturbatively, leaving the Regge slope unchanged and leading to a moderate change of the intercept. The resulting spectra are in very good agreement with experiment [49,76]. The situation is more intricate for the case of heavy quarks, like c-quarks, since now conformal symmetry is strongly broken and the occurrence of linear trajectories is far from obvious. It has been shown, however, that the formalism can indeed be extended to heavy quark bound states [77,78], leading to a fair agreement with the data. In that case, the Regge trajectories are still linear, but the slope depends on the heavy quark mass. The intercept changes quite drastically with the quark mass.
The J/ψ trajectory obtained in [78] is where κ c = √ λ c = 0.874 GeV. This result agrees with the one obtained in a phenomenological potential model [79]. The large change of the intercept as compared to light quarks removes the small-x singularity of quark distribution functions while keeping the counting rules at large Q 2 and at large x unchanged. The change of the slope affects directly only the generalized parton distribution function (GPD). The validity of any mass correction approach should be re-examined carefully in this case. Fortunately, the distribution difference [c(x) −c(x)] is not sensitive to the choice of mass correction methods, since the quark mass affects charm and anticharm distributions in the same way.
In practice, one needs to truncate the expansion in Eq. (7) to have numerical results. For simplicity, we only keep the first term, which is the lowest Fock state containing the charm quark components. The coefficient c τ is determined by the lattice results of G c E (Q 2 ) and G c M (Q 2 ) at the physical limit. We perform a fit to the extracted results of G c E (Q 2 )| physical and G c M (Q 2 )| physical , i.e., the bands in Figs. 1 and 2. As the lattice data from different ensembles are evaluated at different Q 2 values, and have been utilized to determine the quark mass, lattice spacing, and finite volume effects, the effective number of data points in the physical limit is 6 for G c E (Q 2 )| physical and 6 for G c M (Q 2 )| physical 1 . To really capture the uncertainty, we create 200 replicas from the extracted bands. Each replica is firstly generated by randomly sampling 6 data points of G c E (Q 2 )| physical and 6 data points of G c M (Q 2 )| physical from the extracted bands within 0 < Q 2 < 1.4 GeV 2 , which are covered by the lattice data. Then for each data point, the central value is resampled with a Gaussian distribution according to its uncertainty. In addition, we also randomly shift the value of κ c within ±5% in each single fit to one replica to incorporate the theoretical uncertainty. The coefficient determined from the fit is c 5 = 0.018 (3).
We use the same universal form of the function w(x) from [54] with a = 0.480 [80]: It incorporates Regge behavior at small x with the J/ψ intercept (13), w(x) → x as x → 0, and the inclusive-exclusive counting rule at large x, q τ (x) → (1 − x) 2τ−3 , as x → 1 [54]. We obtain the asymmetric charm-anticharm distribution function Fig. 3. The result from the fit is in agreement with the qualitative analysis at the beginning of this section that the charm quark tends to carry larger momentum than the anticharm quark based on the charm quark form factors from the lattice calculation.  (15). (15) The [c(x) −c(x)] distribution result is about 3 times smaller in magnitude than the s(x)−s(x) distribution obtained with the same formalism [46]. Although a small asymmetry could be a result of the cancellation of two relatively large c(x) andc(x) distributions, it is possible that the intrinsic charm and anticharm distributions are both small. Furthermore, the charm and anticharm distributions at high energy scales are dominated by the extrinsic sea from perturbative radiations. The experimental observation and isolation of the intrinsic charm effect are extremely challenging in such cases. Thus it is not surprising that the recent measurement of J/ψ and D 0 productions by the LHCb collaboration [13] found no intrinsic charm effect. An ideal place to investigate the intrinsic charm would be the J/ψ or open charm productions at relatively low energies, e.g., at JLab, although it is also possible to see intrinsic charm effects in very accurate measurements of high energy reactions. In addition, lepton-nucleon scattering may provide a cleaner probe than nucleon-nucleon scattering to help reduce backgrounds and increase the chance to observe the intrinsic charm effect, and therefore the future EIC will provide such opportunities.
The nonzero value of G c E (Q 2 ) can also originate from the interference of the q → gq → ccq and q → ggq → ccq sub-processes, without the existence of IC. However, as mentioned earlier, this extrinsic [c(x) −c(x)] asymmetry which arises at the next-tonext-to-leading order level is negligible [38]. Moreover, according to [38], this extrinsic asymmetry would result in a much smaller and negative value of the first moment of [c(x) −c(x)] distribution x c−c compared to x c−c = 0.00047 (15) [46], the EMC measurement [8], and perturbative QCD computation [38] seem to indicate extremely small values of extrinsic charm for x > 0.1. The present determination of the [c(x) −c(x)] distribution gives a strong evidence from LQCD for the existence of nonperturbative intrinsic heavy quarks in the nucleon wavefunction at large x ∼ 0.4 − 0.5 with a magnitude consistent with experimental signals. A consequence of this result is Higgs production at large x F > 0.8 in pp collisions at the LHC from the direct coupling of the Higgs to the intrinsic heavy quark pair [81].

Conclusion and outlook
In this article, we have presented the first lattice QCD calculation of the charm quark electromagnetic form factors in the physical limit. This first lattice QCD calculation indicates that a nonzero charm electric form factor corresponds to the intrinsic charmanticharm asymmetry in the nucleon sea, thereby providing an indication of the existence of nonzero intrinsic charm based on a first-principles calculation. In addition, the nonzero value of the charm magnetic form factor indicates a nonzero orbital angular momentum contribution to the nucleon coming from the charm quarks. We have discussed that the existence of IC is supported by QCD and how an accurate knowledge of the intrinsic charm can help to remove bias in the global fits of PDFs and related phenomenological studies.
Motivated by the new lattice results, we have used the nonperturbative light-front holographic framework incorporating the QCD inclusive-exclusive connection at large x to determine the [c(x) −c(x)] asymmetry up to a normalization factor, which is constrained by the lattice QCD calculation. Since the LFHQCD calculation starts from a nucleon Fock state with hidden charm, the parton distributions determined in this model refer exclusively to intrinsic charm where the small-x behavior is determined by the J/ψ intercept. On the other hand, contributions from gluon splitting are supposed to be determined by the pomeron trajectory with a much higher intercept. These features will be discussed in a separate publication.
The new determination of the [c(x) −c(x)] asymmetry presented here gives additional elements and further insights into the existence of intrinsic charm. It also can provide complementary information to the global fits of PDFs which look for the possibility of IC in the absence of ample experimental data.