Parton distribution functions and constraints on the intrinsic charm content of the proton using BHPS approach

In this work, a new set of parton distribution functions taking into account the intrinsic charm (IC) contribution is presented. We focus on the impact of the EMC measurements on the large $x$ charm structure function as the strongest evidence for the intrinsic charm when combined with the HERA, SLAC and BCDMS data. The main goal of this paper is the simultaneous determination of the intrinsic charm probability ${P}_{c{\bar c/p}}$ and strong coupling $\alpha_s$. This allows us to study the interaction of these two quantities as well as the influence on the PDFs in the presence of IC contributions. By considering $\alpha_s$ which can be fix or free parameter from our QCD analysis, we find that although there is not a significant change in the extracted central value of PDFs and their uncertainties, the obtained value of ${P}_{c{\bar c/p}}$ change by factor 7.8\%. The extracted value of ${P}_{c{\bar c/p}}$ in the present QCD analysis is consistent with the recent reported upper limit of $1.93\%$, which is obtained for the first time from LHC measurements. We show the intrinsic charm probability is sensitive to the strong coupling constant and also the charm mass. The extracted value of the strong coupling constant $\alpha_s(M_Z^2) = 0.1191 \pm 0.0008$ at NLO is in good agreement with world average value and available theoretical models.


I. INTRODUCTION
The distribution of heavy quarks in the proton can provide a comprehensive overview of the nucleon structure, and it is essential for processes in the perturbative QCD (pQCD) calculations at the LHC physics.
According to the Brodsky, Hoyer, Peterson, and Sakai (BHPS) model [1,2], the charm quark distribution of the nucleon consists of two separate contributions. The first one is generated perturbatively by gluon splitting to cc (g → cc) from DGLAP evaluation equations [3], * Abdolmalki@semnan.ac.ir † Khorramiana@semnan.ac.ir which is generally referred to "extrinsic charm", see [4][5][6]. This extrinsic charm distribution depends logarithmically on the charm mass m c and is most important at low x. Next is the "intrinsic charm" (IC) which has a nonperturbative origin, in contrast to extrinsic charm, and is associated with bound state hadron dynamics. This distribution is described by Fock state in the nucleon structure which is dominant at large x arises from multiply connected charm quark pairs by the valence quarks.
With same rapidities of the quarks which constituent the nucleons, this contribution is maximal in the minimal offshellness and depends on the non-perturbative structure of the nucleon. While the intrinsic charm contribution is dominant at high x and depends on 1/m 2 c , the extrinsic charm distribution is dominant at small x and depends logarithmically on charm mass.
The intrinsic charm contribution is important for estimating the flux of high energy neutrino which observed in the IceCube measurement. In Ref. [13], the prompt neutrino spectrum using the intrinsic charm contribution is estimated. It is shown that the prompt atmospheric neutrino flux taking into account the intrinsic charm contribution is comparable with the extracted results of QCD calculations without taking to account the intrinsic charm. Undoubtedly, the IceCube measurements will constrain the IC component in the proton. This kind of measurements will also contribute to the main questions in high energy physics phenomenology in the future [42].
Not only intrinsic charm quark but also intrinsic strange and bottom quarks are principle property of the wave functions of hadronic bound states. The application of the operator product expansion shows that the probability for the analogous intrinsic bottom contribution from |uudbb > Fock states is suppressed relative to intrinsic charm by a factor of m 2 c /m 2 b . There are some processes which are sensitive to the charm quark distribution in the large x region. For example the European Muon Collaboration (EMC) [43] provided the charm structure function data F c 2 as strong evidence for intrinsic component at large x [43]. The measurement at x=0.42, Q 2 = 75 GeV 2 is approximately 30 times higher than predicted from gluon splitting. It should be noted that the EMC data is the only DIS evidence for intrinsic charm at high x and it is worthwhile to include these data in the study of intrinsic charm.
A description of the EMC data with intrinsic charm is presented in Refs. [38,39,44,45]. In Ref. [38], the NNPDF included the EMC data with using the standard W 2 > 12.5 GeV 2 cut to avoid the presence of dynamical higher twists. Although the EMC data are included without nuclear corrections in NNPDF collaboration [38], several fit variants have been performed to assess the impact of nuclear corrections on the EMC data and find a moderate impact.
In Ref. [44], Hoffmann and Moore performed the NLO calculation of the BHPS model using EMC data. The first NLO analysis of both extrinsic and intrinsic contributions with EMC data was done by Harris, Smith, and Vogt [45].
Since the first experimental evidence of IC originated from the EMC measurement at large x, a variety of charm hadrons measurements are consistent with the existence of IC. In the case of hadron-hadron collisions, the intrinsic charm leads to the production of charm hadrons such as the Λ c (cud) as observed in ISRF experiments [46] and more recently by SELEX [47] at high x F = x c + x u + x d from the coalescence of the charm quark with its coming valence quarks, as well as quarkonium productions at high x F = x c + xc. Double intrinsic charm Fock states such as |uudcccc > in the proton lead to the production of double-charm baryons as well as double quarkonia at high x F , see Ref. [48]. Previous fixed target J/ψ measurements also show signs of significant IC contribution, taking into account the nuclear mass dependence, as measured at CERN and Fermilab, see e.g. [49,50].
On the other hand, the production of the prompt photons at hadron colliders in pp(p) → γ + c-jet is sensitive to the charm distribution and may provide good evidence for IC at high p T . In Refs. [51][52][53] the results of c-jet production accompanied by vector bosons Z, W ± is stud-ied by using the intrinsic charm quark component. The impact of IC contributions for the prediction of γ + cjet production in pp collisions at the LHC is reported in Ref. [54], where this impact for large p γ T can be distinguished from the case in which the IC contribution is not considered. Without considering IC contribution, the D0 data for γ + c-jet differential cross section [55] can not be described by solely using an extrinsic charm PDF based on the DGLAP evolution equation. In contrast, the D0 data for γ + b-jet differential cross section is explained by the extrinsic standard PDFs. This is consistent with the fact that the ratio of the intrinsic bottom distribution to the intrinsic charm distribution in the proton is suppressed by m 2 c /m 2 b . In addition to the typical observables for IC, the intrinsic bottom (IB) quarks also contribute to diffractive Higgs production in which the Higgs boson carries a remarkable fraction of the proton momentum. Measuring the high x F process pp → HX via intrinsic quarks (IQ) at the LHC would give new constraints on the Higgs couplings to quark pairs, including H → bb [56].
Additionally, the new ATLAS measurements at the LHC [57] have shown that the production of prompt photons accompanied by a charm-quark jet in pp collisions are sensitive to the intrinsic charm content of the nucleon. By having these new ATLAS measurements at 8 TeV, the upper limit of the intrinsic charm probability is found to be about P cc/p =1.93% [58]. They proposed a method that reduced the uncertainty on the determination of intrinsic charm probability. Undoubtedly, demonstrating the compatibility of the DIS data and the new LHC data for the P cc/p extraction [58], would be worth. Regarding the importance of intrinsic heavy quarks, now we have sufficient motivation to incorporate the IC contributions in our QCD analysis.
In this paper, we perform our QCD analysis based on the xFitter open source framework [59,60], which was previously known as HERAfitter. Recently, the low x resummation of QCD analysis is performed by xFitter [61], which leads to the better description of the data at low x and low Q 2 . Other QCD analysis based on xFitter framework are performed in Ref. [62,63]. In Refs. [62,64] we extracted the strong coupling constant using HERA I and II combined data and Neutrino-nucleon structure function data using xFitter. More recently, we also determined the strong coupling constant with polarized data [65]. Note that one of the main purposes of this paper is to determine α s concurrently with the IC probability.
The paper is organized as follows. In the Sec.II, we briefly outline the basic formalism and provide a theoretical overview of the intrinsic heavy quark distributions. In this section, we also introduce the Q 2 evolution of the intrinsic heavy quark distribution functions and the intrinsic heavy quark component of the structure functions. The experimental data which we use in the present analysis is presented in Sec.III. The PDF parametrizations are discussed in Sec.IV. In Sec.V, we present the fit results for the PDFs including the intrinsic charm at NLO, and we make comparisons with other different theoretical analyses from literature. Finally, our discussion and conclusion are given in Sec.VI.

II. BHPS MODEL AND EVOLUTION OF INTRINSIC HEAVY QUARK STRUCTURE FUNCTION
According to the light-cone formalism, the proton wave function can be represented as a sum over the complete basis of free quark and gluon states based on a superposition of the proton wave function and n-particle Fock fluctuation components, |q 1 q 2 q 3 , |q 1 q 2 q 3 g , |q 1 q 2 q 3 qq , etc, [1] as where the light-front wave functions ψ j/p (x i , k ⊥i ) depends on the relative momentum coordinates k ⊥i and x i = k + i /P + in which k i explains the parton momenta and P denotes the hadron momentum. For a proton with mass of m p , the general form of the Fock state wave function is where, m 2 i = m 2 i + k 2 ⊥i is the square of the average transverse mass of parton i. Note that by decreasing m 2 p − M 2 , the Γ as a vertex function, expected to be a slowly varying. This form can be applied to the higher Fock components for an arbitrary number of light and heavy quarks. The momentum conservation for n number of partons in state |j demands n i=1 k ⊥i = 0 and n i=1 x i = 1. In the BHPS model [1], the probability distribution in the 5-particle intrinsic Fock state is where N 5 normalizes the 5-particle Fock state probability that is determined by P cc/p =´1 0 dx 1 ...dx 5 dPIC dx1,...,dx5 , where P cc/p is the intrinsic charm probability in the proton. In the heavy quark limit, m c , m c ≫ m p , m q , where the mass of light quark and proton are negligible compared to the heavy quarks mass, the Fock state probability distribution is written as where N 5 = N/m 4 c,c is determined from Eq. (4) by integrating over dx 1 ...dxc, so the intrinsic heavy quark prob-ability distribution is given by We can simplify the above equation by considering P cc/p = 1% IC probability. For more details see Refs. [23][24][25].
According to Eq. (5), as the BHPS model predicts [1], the existence of the intrinsic bottom (IB) distribution is very similar to the IC, but differs in the normalization factor by a coefficient m 2 c /m 2 b . Therefore, P bb/p = P cc/p (m 2 c /m 2 b ) ∼ 0.001 taking into account the 1% IC probability distribution.
As mentioned above, the heavy quark distribution in the standard approach, which is used by almost all global PDFs analyses, is generated by quark-gluon fusion in the DGLAP equations at the starting scale on the order of the heavy quark mass. As the BHPS model predicts, the purely perturbative treatment can not give a good description of the proton structure. Therefore the full charm parton distribution must be expressed by the sum xc(x, Q 2 ) = xc ext (x, Q 2 ) + xc int (x, Q 2 ), where xc ext (x, Q 2 ) and xc int (x, Q 2 ) indicate the extrinsic charm that is radiatively generate by the DGLAP equation (perturbative) and the intrinsic charm (non-perturbative) at an initial scale Q 0 ≃ m c , respectively. This decomposition is a good approximation at any scale since the intrinsic charm is controlled by non-singlet evaluation equations [66]. Therefore the evaluation of heavy quarks must be divided into independent parts, singlet and nonsinglet. The compact form of the DGLAP equation used by the standard global analyses, is given bẏ whereḟ i denotes the light quarks, heavy quarks and gluon and P ij is the splitting function which is known up to three-loop order in the massless M S scheme [67,68]. By considering the intrinsic heavy quark distribution Q int in the proton, one can find the non-singlet equation as [66]Q In perturbative QCD, the structure functions can be written as a convolution between the hard scattering coefficient function and parton distribution functions which are parametrised and determined from experimental data.
The deep-inelastic structure functions F (x, Q 2 ) consists of the heavy and light flavor component where F h is the heavy contribution of DIS structure function, which is valid if only the heavy quark electric charge is non-zero.
As we plan to study the intrinsic charm component of the proton structure function, then we need to add the intrinsic structure function F (int) (x, Q 2 ) to the extrinsic heavy structure function F h,RT (ext) (x, Q 2 ) using the Thorne-Roberts (RT) scheme [69] as The definition of the heavy quark probability distribution in Eq. (5), using P cc/p = 1%, for the |uudcc Fock state, leads to a simple form for the IC component in the proton structure function, as F cc/p 2(int) (x, Q 2 ) = 8 9 x´dx 1 · · · dxc dPIC dxi···dxcdxc = 8/9xc int (x, Q 2 ) with c int from Eq. (5), when the charmed mass is negligible in the leading order [7]. By considering mass effects, the IC component in the proton structure F cc/p with c(z, γ) = c int (z) − zc int (γ)/γ for z ≤ γ. Also, γ as another mass scaling variable which depends on Q 2 and is defined by In the above,ĝ(ξ, γ) iŝ For more detail see Refs. [44,70].

III. DATA SETS
For the QCD analysis of PDFs including intrinsic charm, different DIS processes which are generally important in the presence of intrinsic charm are used to determine the unknown parameters in the PDFs parametrization, strong coupling constant α s and P cc/p as the extra fit parameters.
The neutral current (NC) and charged current (CC) cross section data at HERA have explored a wide kinematic region of the Bjorken variable x and negative boson transverse momentum squared Q 2 . In this paper we used the NC and CC inclusive DIS experimental data collected by H1 ans ZEUS [75][76][77][78][79][80][81][82] from both HERAI and HERAII with the proton beam energy E p = 460, 575, 820 and 920 GeV 2 related to center of mass energy √ s = 225, 251, 300 and 320 GeV, respectively.
High Q 2 CC data, together with difference between the NC e − p and e + p at high Q 2 , constrain the u and dvalence PDFs which dominate at large x [71]. The wide kinematic region of the precise NC and CC cross-sections data allow us to extract PDF sets.
We also choose the proton structure function data which are reported by BCDMS [72] in deep-inelastic scattering of muons on the hydrogen target at the beam energy of 100, 120, 200 and 280 GeV.
The SLAC data on ep scattering used in our analysis, as well. It is obvious that our motivation to choose this data set of experimental data is due to the kinematic range of large x, where the intrinsic charm is dominant.
In the HERA kinematic domain, where the virtuality Q 2 of the exchanged boson is small, Q 2 ≪ M 2 Z , charm production is dominated by virtual photon exchange, where the xF cc 3 (x, Q 2 ) is negligible. Therefore the neutral current deep-inelastic ep cross section by considering the running electromagnetic coupling, α(Q 2 ), without QED and electro-weak radiative corrections, may be written in terms of the structure functions F cc 2 (x, Q 2 ) and F cc L (x, Q 2 ). Since the F cc L contribution is suppressed only at low rapidity y [86] at HERA, it is possible to neglect this contribution. So the important contribution to the reduced charm cross section is the charm structure function, F cc 2 . Therefore additional correction must be done in this contribution in the presence of IC contribution. The importance of this kind of data set is due to the existence of charm PDFs, where the IC contribution can be included.
Since the main goal of this paper is to determine the intrinsic charm probability in the proton, we need to include the charm structure function data in the high x region, where the intrinsic charm is dominant. Therefore, in this analysis, we choose the EMC measurement of the large x charm structure function as a first experimental evidence of IC. These data are produced in inclusive dimuon and trimuon iron target. Since the EMC F c 2 data are extracted using the nuclear target, we need to consider the nuclear corrections in the present QCD analysis, as we will explain in the next section.
Although the fixed-target Drell-Yan data [87] are sensitive to the distributions of anti-quark in large x and can be used in the QCD global analysis, we did not find a significant impact on PDF parametrization by including this kind of dataset in the presence of IC contribution. On the other hand, in our previous analysis [63], we have also investigated that including the jet and W , Z boson data at the LHC [88][89][90][91] do not any impacts on the PDFs in the presence of IC contribution. So we can exclude these data sets in the present analysis.
In this analysis, we apply the (W 2 − W 2 th )/W 2 with W 2 th = 16 on heavy quark structure function to suppress charm production near threshold as suggested in Ref. [92]. Note that, we exclude the data using Q 2 < 4 GeV 2 and W < 15 GeV 2 cuts, where the higher twist correction might become relevant. So the total number of data points in our QCD analysis with considering these cuts are reduced from 1939 to 1535.

IV. PDF PARAMETRIZATION
In this article, we perform our QCD analysis using the xFitter open source framework [59,60], which previously was known as HERAfitter [93]. Very recent analyses using xFitter package are reported in Refs. [61,62].
The numerical solutions of the DGLAP evolution equations for PDFs in pQCD framework at NLO are implemented in the QCDNUM package in x-space [94]. To compute the perturbative part of heavy quark contributions of the DIS structure-function, we used the optimal Thorne-Roberts (RT-opt) [69], General-Mass Variable Flavor Number (GM-VFN) scheme. This scheme assumes that the charm distribution is generated perturbatively by gluon and light quark splittings, and its value depends strongly on the charm mass. If the PDFs are parameterized as a function of x at initial scale Q 2 0 , the QCD evolution will help us to achieve it at any value of Q 2 . In the present analysis, the initial QCD scale is chosen to be Q 2 0 = 1.9 GeV 2 , which is below the charm threshold. In this approach, we choose the heavy quark masses m c = 1.50 GeV and m b = 4.5 GeV, and we take µ r = µ f = Q for the QCD re-normalization and factorization scale.
We choose a standard parametrization for the PDFs at the input scale Q 2 0 = 1.9 GeV 2 [71] to be: The parametrised PDFs are the valence distributions, xu v , xd v , the u-type and d-type anti-quarks distributions, xŪ , xD, and the gluon distribution, xg. We assume the relations xŪ = xū and xD = xd + xs at the starting scale.
In summary, our parametrisation is: In the above equations, to ensure the same behavior of the xū and xd as x → 0, one can impose the additional constraints BŪ = BD and AŪ = AD(1 − f s ) [71]. In this analysis, we fixed the DŪ parameter after the first minimization because, in the presence of IC contribution, the selected DIS data set does not constrain the DŪ well enough. The strange-quark distribution is expressed as an x-independent fraction, f s , of the d-type sea, xs = f s xD, at Q 2 0 . The value f s = 0.31 was chosen in the strange quark density as suggested in Ref. [95]. For the gluon PDF, C ′ g is fixed to C ′ g = 25 to ensure a positive gluon density at large x, as suggested in Ref. [96].
The normalization parameters for gluon and valence distributions, A g , A uv and A dv are constrained by the fermion number and momentum sum rules, In fact, in the above sum rule the total intrinsic charm quark momentum fraction iŝ Notice that in Eq. (5) we can fix the P cc/p value with 1%, which is only included in the intrinsic part of above equation, or can be considered as a free parameter. It should be noted that the extrinsic charm generated perturbatively using the DGLAP evaluation equation and the intrinsic charm evolve by a non-singlet evaluation equation [54].
In this analysis, we use the nuclear correction on the EMC data because these data come from a nuclear target. This correction creates a connection between the parton distributions in the nucleus A and parton distribution in the proton which we model as: where is the parton distribution with type i in the nucleus and f i (x, Q 2 ) is the corresponding parton distribution in the proton. A and Z are the mass number and atomic number, respectively. To study the impact of the EMC data on the PDFs behavior considering IC contributions, we apply the nuclear correction factor from Ref. [97] on EMC data.
By comparing the theoretical and experimental measurements of various physical observables, we determine the unknown PDFs parameters by minimizing the χ 2 function taking into account the correlated and uncorrelated measurement uncertainties. The χ 2 function is minimized using the MINUIT package [98], and is defined as [61,99] where C −1 ij is the covariance matrix. If the full correlated uncertainties of the experimental data are available, the χ 2 function is as follows where t i is the theoretical prediction, d i is the measured value of the i-th data point, (Q 2 , x, s), and δ i,stat and δ i,unc are the relative statistical and uncorrelated systematic uncertainties. In the above, β i j are the corresponding systematic uncertainties, in which case s j are the nuisance parameters associated with the correlated systematic error. Here j labels the sources of correlated systematic uncertainties, and in the Hessian method the s j are not fixed. If the s j are fixed to zero, the correlated systematic errors are ignored.
The Hessian method for the PDF uncertainties are obtained from ∆χ 2 = T 2 . A tolerance parameter T is selected such that the criterion ∆χ 2 = T 2 ensures that each data set is described within the desired confidence level. The correlated statistical error on any given quantity q is then obtained from the standard error propagation: The Hessian matrix is defined as H α,β = 1 2 ∂ 2 χ 2 /∂p α ∂p β , and thus the covariance matrix C = H −1 is the inverse of the Hessian matrix evaluated at the χ 2 minimum. In order to be able to calculate the fully correlated 1σ error bands for the valence PDFs, one can choose T = 1 in the xFitter package.
Indeed, we have 13 unknown PDFs parameters in addition to free parameters α s and P cc/p which are obtained from the fit.

V. FIT RESULTS
The main goal of this paper is the simultaneous determination of the intrinsic charm probability P cc/p and strong coupling α s . This allows us to study the interaction of these two quantities as well as the influence on the PDFs in the presence of IC contributions.
For fit analysis, we will include the DIS reduced cross section and differential cross sections data from HERA I+II, the charm combined cross sections data from H1 and ZEUS, and the proton structure function data from BCDMS, SLAC, and EMC. The detailed information and references, number of data points, and their kinematic range of x and Q 2 for each data set are summarized in Table I.
To investigate the effect of the EMC data in the present analysis with and without IC contributions, we will divide our QCD analysis into four fits: • Base: We include all the data sets of Table I with exception of the EMC experimental data, this totals 1519 data points. In the Base fit we consider zero IC contribution, i. e., P cc/p ≡ 0. We fixed α s (M 2 Z ) = 0.1182 considering to the world average data as a first step.
• Fit A: We include all the data from the Base fit in addition to the EMC F c 2 data and this gives us 1535 data point. This fit also assumes zero IC contribution and we use a fixed α s (M 2 Z ) = 0.1182.
• Fit B: We use all the data sets of the Table I and we use a fix α s (M 2 Z ) = 0.1182 parameter and include an intrinsic charm probability P cc/p value as a free parameter.
• Fit C: This is all same as Fit B, except now use a free α s (M 2 Z ).
The χ 2 /number for each data set after cuts, the total χ 2 , and χ 2 /dof for all fits are summarized in Table II. According to Table II and without taking into account IC contribution in Base fit, as a first step, the total χ 2 /dof value is 1872/1506 = 1.243. There are 13 unknown parameters for PDFs only.
As a second step and the same as Base fit, we have 13 unknown parameters for PDFs in Fit. A. In this case and according to Table II, we obtain 2047/1522 = 1.345 for the total χ 2 /dof for Fit A.
For a more complete discussion, the investigation of necessity to include a non-zero P cc/p in the present QCD analysis and also the effect of including the EMC data would be important. In Fit B, as a third step, we obtain 1959/1521 = 1.289 for the total of χ 2 /dof for Fit B. There are 13 unknown parameters for PDFs and an unknown parameter for intrinsic charm probability, therefore the total number of unknown parameters is 14.
Our motivation to present Base, Fit A and Fit B results is to show that an IC component is unnecessary as long as the EMC data remains excluded.
Finally, in Fit C we consider the intrinsic charm probability P cc/p and α s (M 2 Z ) as free parameters in our QCD analysis, so there are 13 unknown parameters for PDFs and two unknown parameters for strong coupling constant and intrinsic charm probability. In Fit C, the total number of unknown parameters is 15, which should be obtained by the QCD fit on experimental data. In this fit, we obtain 1957/1520 = 1.287 for the total χ 2 /dof with considering intrinsic charm probability P cc/p and α s (M 2 Z ), as the free parameters. It should be noted that, in the present analysis and for Fit A, B, and C, the nuclear effects are considered.
Although the Fit C contains the complete analysis, the comparison of this fit with other above cases with each other would be interested.
Obviously, in comparison Fit A with Fits B and C, there are almost 4% improvements for χ 2 values and the fit quality.  In Ref. [27,100], the NNPDF collaboration studied the influence of the EMC data in their analyses with and without the IC contribution. They found χ 2 per point = 7.3 for the EMC data if the charm PDF was generated perturbatively, and χ 2 per point = 4.8 when the IC contribution was included and the EMC data was rescaled. They improved the fit quality of the EMC data by imposing an additional cut of x > 0.1 when the charm PDF was fitted.
According to Table. II, we note that the χ 2 per point for the EMC data are 6.37, 4.56 and 4.37 in Fits A, B and C, respectively. When the IC includes in Fits B and C, we find an improvement in the χ 2 per point for the EMC data of 28% and 31%, respectively in comparison with Fit A.
To investigate the specific impact of the EMC data with and without IC contribution at NLO, we need to compare our results for our individual fits: Base, Fits A, B, and C. In Fig. 1, we display some samples of our theoretical predictions for all fits in comparison to fixed target DIS data from HERAI+II [71], H1-ZEUS combined charm cross section [74] (Fig.1-a), BCDMS [72], SLAC [73] and EMC [43] measurements ( Fig.1-b) and their uncertainties at NLO as a function of x for different values of Q 2 . According to Fig. 1(b), the data description turns out to be better including the IC component in comparison with the Fit. A results. In fact, we find that in general, it is necessary to include a non-zero P cc/p to have a better description of the EMC data.
In Table III, the numerical QCD fit results for the PDFs parametrization according to Eq. (13), α s value and the intrinsic charm probability P cc/p are summarized in four separate cases.
It is obvious that the P cc/p value has significant sensitivity to the α s value. The simultaneous determination of P cc/p and α s , as shown in Table III, gives a somewhat lower IC probability and an α s which is more or less in agreement with the world average value of α s (M 2 Z ) = 0.1182 ± 0.0011 [101]. This is a very reasonable result based on the fact that IC is a non-perturbative phenomenon.
According to Tables II and III, the increase of 8% of χ 2 /dof in Fit A with respect to our Base fit, and ∼ 4% improvement of χ 2 /dof in Fit B and Fit C with respect to Fit A, is obtained. Conversely, we observed the significant changes in both, the PDFs parameters and their uncertainties. Using the definition of ∆P cc/p = P F it.C cc/p − P F it.B cc/p , we find a change of ∼7.6% on ∆P cc/p /P F it.B cc/p . It is clear that a simultaneous determination of P cc/p and α s in the present analysis can impact the central value of the intrinsic charm probability. Fig. 2 illustrates the NLO QCD fit results for valence, sea and gluon PDFs as a function of x at the initial scale of Q 2 0 = 1.9 GeV 2 . Although there are no significant changes in the valence and sea PDFs in all fits, significant changes in the gluon PDFs are observed (left panels). To clarify the difference between our four fits, we present the relative uncertainties δxq(x, Q 2 )/xq(x, Q 2 ) with respect to the Base fit. According to this figure, the changes of the central values of the PDFs, their uncertainties, or both are observed (right panels).
In Fig. 3, we present our results for the extrinsic charm PDF based on our four fits as a function of x for Q 2 = 10 GeV 2 and Q 2 = 100 GeV 2 at NLO (left panels). The relative uncertainties δxc(x, Q 2 )/xc(x, Q 2 ) are also shown (right panels).
Since our analysis fixed the charm mass to 1.5 GeV, we are interested in study the impact of the charm mass value on P cc/p . Basically, the charm mass impact in our analysis has led us to recalculate the simultaneous fit of P cc/p and α s for different values of charm mass. Here, we choose Fit C, as a reference point to study the impact of charm mass on PDFs behavior and the extraction P cc/p and α s values.
As mentioned above, to compute the perturbative part of structure-function, we utilize the GM-VFN scheme, which assumes that the charm distribution is generated perturbatively by gluon and light quark splittings. The detail of this is sensitive to the charm quark mass, which is not precisely known. We are going to use the RT-opt heavy quark scheme and we are going to vary the charm mass 1.2 GeV to 1.8 GeV and we will scan the χ 2 to determine the optimal value [74]. In the present paper, we consider the initial value of Q 2 0 for the evolution to be below the charm threshold µ c .
The χ 2 (m c ) value is calculated for each fits and the optimal values of m opt c parameter is determined in a given scheme from a parabolic fit to the χ 2 (m c ) values as where χ 2 min is the χ 2 value at the minimum and σ(m opt c ) is the fitted experimental uncertainty on m opt c . The procedure of the χ 2 -scan is illustrated in Fig. 4 for the optimal RT (RT-opt) scheme by fitting the experimental data. According to this figure, we find the minimum of χ 2 value for the charm mass of 1.50 GeV. This was our motivation to choose m c = 1.50 GeV in the presence of IC.
In Table IV, we summarize the influence of the charm mass for the determination of P cc/p and the strong cou-pling constant using Fit B and Fit C. According to this table and for Fit B, P cc/p goes from 0.78% to 1.4% as increases m c from 1.39 to 1.6 GeV. For Fit C, we found P cc/p changes from 0.97% to 0.88% as m c changes and we also find the α s (M 2 Z ) increase from 0.1152 to 0.1225 with the similar uncertainties.
As expected, we show the intrinsic charm probability and the strong coupling constant values depend on the charm mass value. The stronger dependence of P cc/p and charm mass is due to the fact that the heavier m c makes it harder to create cc pairs to obtain 5-particle Fock states.
In Since the gluon PDF plays an important role in the total cross section of the top quark production, we used the results of Fit C with different charm quark masses and Fit B with fix m c = 1.5 to calculate the top quark cross section at the LHC using the HATHOR package [104].
According to Fig. 5, the gluon densities show significant variation for the different cases with varying charm mass values for Fit C, as a complete fit to study the impact of charm mass. We obtain the total top quark cross section as 804.7±6. 6 Fig. 6 shows the NLO extrinsic charm PDF extracted from Fit C with different charm quark mass values for range of Q 2 as a function of x, for different values of Q 2 = 10, 100 GeV 2 (left panels). The relative uncertainty ratios xc(x, Q 2 )/xc(x, Q 2 ) ref in respect to m c =1.5 GeV (right panels) are shown.
In Table V, we present the intrinsic charm probability values extracted from Fit C with the present analysis for different values of charm mass and compare with other predictions from the literature. The extracted value of intrinsic charm probability is in good agreement with the other predictions of the models.
In Table VI, we present the momentum fraction of the IC distribution < x > int c+c according to Eq. (16) and total charm momentum fraction, < x > tot c+c =< x > int c+c + < x > ext c+c , extracted in Fit C, and compared to the NNPDF3IC and CT14-IC at Q 2 = 4 GeV 2 . Our results  for < x > int c+c and < x > tot c+c are obtained using P cc/p (%)= 0.94 ± 0.20. The momentum fraction of the IC distribution and the total momentum fraction of the charm PDF are in good agreement with other models. Note that in the CT analysis several models for IC are studied whilst NNPDF3IC assumed a specific parametrisation of the charm PDF which is fitted simultaneously to light quarks and gluons.
In Fig. 7, we show the intrinsic charm distribution xc int (x, Q 2 ) with its uncertainty as a function of x at different values of Q 2 with P cc/p = 0.94±0.2%. In Fig. 8, we display the intrinsic charm xc int , extrinsic charm xc ext and the total charm xc int + xc ext distributions with their uncertainties with P cc/p = 0.94 ± 0.2%, as a function of x and Q 2 .
In Fig. 9, we display the valence, sea and gluon PDFs extracted in our QCD Fit C for selected Q 2 values as a function of x, and compared them with the results obtained by HERAPDF2 [71] and NNPDF3IC [38]. Note that the NNPDF collaboration used a different methodology to parametrize the PDFs, and chose different input scales and the cut values.

VI. DISCUSSION AND CONCLUSION
We performed a QCD analysis using fixed target DIS data from HERAI+II [71], H1-ZEUS combined charm cross section [74], BCDMS [72], SLAC [73] and EMC [43] considering IC contribution at NLO. We determine the PDFs with their corresponding uncertainties and also the intrinsic charm probability and strong coupling constant using the xFitter package [59,60].
To investigate the effect of the EMC data in the present analysis and to clarify the results, we divided our QCD analysis into four steps. As a first step, we performed a fit with all data sets of Table I, with the exception of the EMC experimental data. We used a fixed α s (M 2 Z ) value of 0.1182 taken from the world average [101] for the Base fit.
We then include all the data from the Base fit in addition to the EMC F c 2 data with assuming zero IC contribution in Fit A.
Taking into account the IC contribution and considering P cc/p as a free parameter, we performed our QCD analysis for two separate Fit B and Fit C. In Fit B, we used all data sets which used in Fit A and we fixed strong coupling α s (M 2 Z ). In this case we consider the intrinsic charm probability P cc/p as a free parameter. We found that P cc/p depends on the α s value, so the simultaneous determination of these parameters in the QCD analysis, taking into account IC is important. In Fit C, we considered P cc/p and α s (M 2 Z ) as free parameters. The extracted value of the strong coupling constant α s (M 2 Z ) = 0.1191 ± 0.0008 at NLO in Fit C, is in good agreement with the world average value of α s (M 2 Z ) = 0.1182 ± 0.0011 [101]. This extracted value is also consistent with our recent reported results of α s (M 2 Z )= 0.1199 ± 0.0031, based on charged current neutrino-nucleon DIS data at NLO [64]. For P cc/p , we extracted 1.017 ± 0.2 for Fit B and 0.94 ± 0.2 for Fit C. We found 7.6% difference between Fit B and Fit C on P cc/p values by comparing these fits.
The Q 2 -evolution of intrinsic charm PDF using P cc/p = 0.94 ± 0.2% in the proton which extracted in the present analysis, are presented. The extracted value of P cc/p in the present QCD analysis is consistent with the recent reported upper limit of 1.93% of Ref. [58], which is obtained for the first time from LHC measurements [57].
To find the impact of the charm mass in the present analysis, we study the influence on the extraction of P cc/p . For Fit B, P cc/p goes from 0.78% to 1.4% as in-creases m c from 1.39 to 1.6 GeV. For Fit C, we also find P cc/p changes from 0.97% to 0.88% as m c changes and the α s (M 2 Z ) increase from 0.1152 to 0.1225. Since the gluon PDF plays an important role in the total cross section of top quark production, we find a total top quark cross section of 817.8±6.4 using by Fit B when m c = 1.5. In Fit C, we find the cross section of 804.7±6.6, 828.7±8.0 and 855.7±7.8 for m c = 1.43, m c = 1.5 and m c = 1.55, respectively.
The parabolic fit to χ 2 as a function of m c is used to find the optimal value of m c . In this analysis, we choose the value of m c equal to 1.50 GeV. We performed our QCD analysis based on this value of m c in presence of IC.
In the comparison of our PDFs to others in the literature, the valence, sea, gluon and charm PDFs and their uncertainties are in good agreement with other theoretical models. This agreement is despite of different phenomenology and differences in the parametrization and choose of the data sets.
In our previous results, we have shown the differential cross section of γ + c-jet is sensitive to intrinsic charm contribution [54]. In the near future, we will investigate the role of IC contributions in other processes, such as γ + c/b, and Z/W + c/b production cross section measurements in pp collisions at the LHC.

ACKNOWLEDGMENTS
We thank S. Brodsky for valuable and constructive suggestions during the planning and development of this research work. We would also like to thank F. Olness, A. Glazov, O. Zenaiev, Z. Karamloo and M. Azizi for helpful discussions, comments and carefully reading the manuscript. A. K. thanks SITP (Stanford Institute for Theoretical Physics) and the Physics Department of SMU (Southern Methodist University) for their hospitality at the beginning of this work. A. K. is also grateful to CERN TH-PH division for the hospitality where a portion of this work was performed. H. A. is also grateful to the DESY FH group for the hospitality and financial support at DESY.  (Fig.1-a), the proton structure function from SLAC [73] and BCDMS [72], and the charm structure function from EMC data [43] (Fig.1-b).      10. The extrinsic charm PDFs from Fit C (left panels) and the relative uncertainties δxq(x, Q 2 )/xq(x, Q 2 ) (right panels) for the selected scales Q 2 = 4, 10 , 100, 8317 GeV 2 , as a function of x, and compared with HERAPDF2 [71] and NNPDF3IC [38].