Constraints on the intrinsic charm content of the proton from recent ATLAS data

We investigate the possibility to constrain the intrinsic charm probability wcc¯=Pcc¯/p\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$w_{\mathrm{c} {\bar{{\mathrm{c}}}}}= P_{{\mathrm{c}{\bar{\mathrm{c}}}} /\mathrm{p}}$$\end{document} using first ATLAS data on the associated production of prompt photons and charm-quark jets in pp collisions at s=8\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{s} = 8$$\end{document} TeV. The upper limit wcc¯<1.93\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$w_{\mathrm{c} {\bar{{\mathrm{c}}}}}< 1.93$$\end{document} % is obtained at the 68 % confidence level. This constraint is primarily determined from the theoretical scale and systematical experimental uncertainties. Suggestions for reducing these uncertainties are discussed. The implications of intrinsic heavy quarks in the proton for future studies at the LHC are also discussed.


Introduction
One of the fundamental predictions of quantum chromodynamics is the existence of Fock states containing heavy quarks at large light-front (LF) momentum fraction x in the LF wavefunctions of hadrons [1,2]. A key example is the |uudcc intrinsic charm Fock state of the proton's QCD eigenstate generated by cc-pairs which are multiply connected to the valence quarks. The resulting intrinsic charm (IC) distribution c(x, Q 2 ) is maximal at minimal off-shellness; i.e., when all of the quarks in the |uudcc LF Fock state have equal rapidity. Equal rapidity implies that the constituents in the five-quark light-front Fock state have momentum fractions x i = k + i P + ∝ m 2 i + k 2 T i , so that the heavy quarks carry the largest momenta.
The study of the intrinsic heavy quark structure of hadrons provides insight into fundamental aspects of QCD, especially its nonperturbative aspects. The operator product expansion (OPE) predicts that the probability for intrinsic heavy Qquarks in a light hadron scales as κ 2 /M 2 Q due to the twist-6 G 3 μν non-Abelian couplings of QCD [2,3]. Here κ is the a e-mail: gennady.lykasov@cern.ch characteristic mass scale of QCD. In the case of the BHPS model [1] within the MIT bag approach [4], the probability to find a five-quark component |uudcc bound in the nucleon eigenstate is estimated to be in the range 1-2 %. Although there are many phenomenological signals for heavy quarks at high x, the precise value for the intrinsic charm probability w cc = P cc/p in the proton has not as yet been definitively determined. The first evidence for intrinsic charm (IC) in the proton originated from the EMC measurements of the charm structure function c(x, Q 2 ) in deep inelastic muon-proton scattering [5]. The charm distribution measured by the EMC experiment at x bj = 0.42 and Q 2 = 75 GeV was found to be approximately 30 times that expected from the conventional gluon splitting mechanism g → cc. However, as discussed in ref. [6], this signal for (IC) is not conclusive because of the large statistical and systematic uncertainties of the measurement.
A series of experiments at the Intersection Storage Ring (ISR) at CERN, as well as the fixed-target SELEX experiment at Fermilab, measured the production of heavy baryons c and b at high x F in pp, π − p, − p collisions [7][8][9]. For example, the c will be produced at high x F from the excitation of the |uudcc Fock state of the proton and the comoving u, d, and c quarks coalesce. The x F of the produced forward heavy hadron in the pp → c reaction is equal to the sum x F = x u + x d + x c of the light-front momentum fractions of the three quarks in the five-quark Fock state. However, the normalization of the production cross section has sizable uncertainties, and thus one cannot obtain precise quantitative information on the intrinsic charm contribution to the proton charm parton distribution function (PDF) from the ISR and SELEX measurements [10,11]. Another important way to identify IC utilizes the single and double J/ψ hadroproduction at high x F as measured by the NA3 fixed target experiment [12][13][14]. Intrinsic heavy quarks also leads to the production of the Higgs boson at high x F and the LHC energies [15]. It also implies the production of high energy neutrinos from the interactions of high energy cosmic rays in the atmosphere, a reaction which can be measured by the IceCube detector [16].
The first indication for intrinsic charm at a high energy collider was observed in thep p → cγ X reaction at the Tevatron [17][18][19][20][21]. An explanation for the large rate observed for events at high E T based on the intrinsic charm contribution is given in refs. [22,23]. A comprehensive review of the experimental results and global analysis of PDFs with intrinsic charm was presented in [24,25]. Recently, the estimation of w cc = 0.34 ± 0.14% was obtained in the NNPDF3.1 analysis [26] based on the global PDF fitting procedure over a wide range of LHC measurements. In our previous publications [27][28][29][30][31], we showed that the IC signal can be visible in the production of prompt photons and vector bosons Z /W in pp collisions, accompanied by heavy-flavor c/bjets at large transverse momenta and the forward rapidity region (|y| > 1.5), kinematics within the acceptance of the ATLAS and CMS experiments at the LHC. Based on our previous results [27][28][29][30][31], in the present paper we investigate the possibility to constrain the intrinsic charm probability using first ATLAS data [32] on measurements of differential cross sections of isolated prompt photons produced in association with a c-jet in pp collision at √ s = 8 TeV. Note that these ATLAS data were reported recently and were not considered by the NNPDF3.1 group [26]. Thus, the present study could be viewed as an addition to the NNPDF3.1 analysis.

Variable intrinsic charm
Our considerations below are based on two approaches: the first is the analytical QCD calculation and the second is the MC generator Sherpa [33]. We will first present the scheme of our QCD analysis. The systematic uncertainties due to hadronic structure are evaluated using a combined QCD approach, based on the k T -factorization formalism [34,35] in the small-x domain and the assumption of conventional (collinear) QCD factorization at large x. Within this approach, we have employed the k T -factorization formalism to calculate the leading contributions to the O(αα 2 s ) off-shell gluon-gluon fusion g * g * → γ cc. In this way one takes into account the conventional perturbative charm contribution to associated γ c production. In addition there are backgrounds from jet fragmentation. Here α is the electromagnetic coupling constant and α s is the QCD coupling constant.
The IC contribution is computed using the O(αα s ) QCD Compton scattering cg * → γ c amplitude, where the gluons are kept off-shell and incoming quarks are treated as on-shell partons. This is justified by the fact that the IC contribution begins to be visible in the large x (x ≥ 0.1) domain, where its transverse momentum can be safely neglected. The k Tfactorization approach has technical advantages, since one can include higher-order radiative corrections using the transverse momentum dependent (TMD) parton distributions of the proton (see reviews [36] for more information). Technically, we employ the numerical solution of the CCFM gluon evolution equation (see [35] and references therein), which resums the leading logarithmic terms, proportional to log 1/x, up to all orders of perturbation theory. In addition, we take into account several standard pQCD subprocesses involving quarks in the initial state. These are the flavor excitation cq → γ cq, quark-antiquark annihilation qq → γ cc and quark-gluon scattering subprocess qg → γ qcc. These processes become important at large transverse momenta E γ T (or, respectively, at large parton longitudinal momentum fraction x, which is the kinematics needed to produce high E γ T events); it is the domain where the quarks are less suppressed or can even dominate over the gluon density. We rely on the conventional (DGLAP) factorization scheme, which should be reliable in the large-x region. Thus, we apply a combination of two techniques (referred to as a "combined QCD approach") employing each of them in the kinematic regime where it is most suitable. More details can be found in [37] (see also references therein).
According to the BHPS model [1,38], the total charm distribution in a proton is the sum of the extrinsic and intrinsic charm contributions.
The extrinsic quarks and gluons are generated by perturbative QCD on a short-time scale associated within the largetransverse-momentum processes. Their distribution functions satisfy the standard QCD evolution equations. In contrast, the intrinsic quarks and gluons are associated with a bound-state hadron dynamics and thus have a nonperturbative origin. In Eq. 1 the IC weight w cc , i.e., the probability to find the intrinsic cc pair in the proton, is set by the normalization of the xc int (x, μ 2 0 ) distribution: The total distribution xc(x, μ 2 0 ) satisfies the QCD sum rule, which determines its normalization [39] (see also Eq. 6 in [30]). We define the IC probability as the n = 0 moment of the charm PDF at the scale μ 0 = m c , where m c = 1.29 GeV is the c-quark pole mass.
To get the total distribution xc(x, μ 2 ) at any μ 2 we apply the DGLAP evolution [40][41][42] separately for the extrinsic c ext and intrinsic c int part. As it is shown in [29,30], the interference between the two contributions in Eq. 1 can be neglected, since the IC term xc int (x, μ 2 ) is much smaller than the extrinsic contribution generated at x < 0.1. Note that the evolution in μ 2 equation assumes zero quark mass. The kinematical region, where the IC contribution can be visible in the E γ T -spectrum corresponds to large values of μ ≥ 100 GeV [29,30] which are much larger than the cquark pole mass. Therefore, one can neglect the c-quark mass and apply DGLAP evolution for the IC quark distribution c int . As it is known, there are only few PDFs, e.g., CT14nnloIC or CTEQ66c, which include the IC contribution within the BHPS model [1,2] at fixed values w max cc 1 and 3.5%. In order to get the estimation of the upper limit one needs variable IC contribution to E γ T -spectra at any value of w cc . Therefore, one can adopt a simple relation for any w cc ≤ w max cc [29,30], which provides an interpolation between two charm densities at the scale μ 2 , obtained at w cc = 0 and w cc = w max cc . We have performed a three-point interpolation of all parton (quark and gluon) distributions for w cc = 0, 1 and 3.5 %, which correspond to the CTEQ66 (central value), CTEQ66c1 (BHPS1) and CTEQ66c2 (BHPS2), respectively [10]. For the interpolation function we used the quadratic w cc interpolation, see Eqs. 3 and 4 in Appendix. The difference between quadratic interpolation functions in the interval 0 < w cc ≤ 3.5 % is no greater than 0.5% [29]. Given that, we used the quadratic interpolation for all parton flavors at μ 0 and w cc < w max cc to satisfy the quark and gluon sum rules, see [43,44]. Let us stress that at w cc = w max cc the quark sum rule is satisfied automatically in the used PDF because the intrinsic light qq contributions are included [10]. Note that w cc is treated in xc int (x, μ 2 0 ) of Eq. 1 as a parameter which does not depend on μ 2 . Therefore, its value can be determined from a fit to the data.
In our second approach, we have used the MC generator Sherpa [33] (version 2.2.4) with next-to-leading order (NLO) matrix elements generated by OpenLoops [45] (version 1.3.1) to generate samples for the extraction of the w cc from the ATLAS data. The recent versions of Sherpa generator can provide additional weights [46], which are used to reweight our spectra to PDFs with different IC contribution. In order to obtain weights corresponding to any w we fit three PDF weights: W 0 (0% IC), W 1 corresponding to the BHPS1 (the mean value of the cc fraction is x cc 0.6%, which corresponds to the IC probability w cc = 1.14%) and W 2 corresponding to the BHPS2 ( x cc 2.1%, w cc = 3.54%). The fit is done using a quadratic function on event by event basis and the difference between different fitting functions is negligible (linear and cubic). We used CT14nnloIC [6] PDFs included with the help of LHAPDF6 [47]. The process p + p → γ + any jet (up to 3 additional jets) is simulated with the requirement E γ T > 20 GeV and η γ < 2.7. Additional cuts are applied to match the ATLAS event selection [32]. In order to extract the w-value from the data we first Fwd.
(b) w u. l. = 1.93 % Fig. 1 The E γ T -spectrum calculated with the Sherpa MC generator at NLO compared to the ATLAS data [32]. a Top: the spectrum at central rapidity |η γ | < 1.37 and forward 1.56 ≤ |η γ | < 2.37 region without the IC contribution; a middle: the ratio of the MC calculation to the data for the central rapidity region (w c = 0 %); a bottom: the ratio of the MC calculation to the data for the for forward rapidity regions (w c = 0 %). b The same spectra, as in (a), but with the upper limit of IC contribution w u. l. = 1.93 % validate our E γ T -spectrum generated by the Sherpa event generator in the central rapidity region (|η γ | < 1.37) against γ + c-jet ATLAS spectra. One can see that the ATLAS data are satisfactorily described in the central rapidity region using the Sherpa (NLO) calculation without IC, see Fig. 1 (top).

Results and discussion
The results for the E γ T spectra in the central and forward regions obtained by using CT14nnloIC PDFs without the IC contribution are presented in Fig. 1 (top). One can see from Fig. 1 (top) that the difference between the experimental data and MC calculation are within the total experimental uncertainties, i.e., the experimental ones and theoretical uncertainty due to QCD scale dependence. Therefore, we are unable to determine a precise value of the IC probability from recent ATLAS data. However, one can extract an upper limit for the IC contribution from the data. To obtain this limit we employed a simple method of fitting MC prediction containing IC to the ATLAS data. First, we generate 1000 samples of γ + c E γ T spectra with variable IC approximated by the quadratic formula shown in Eq. (4) in the Appendix. We calculate χ 2 in the forward region for bins 2-11 (counted from zero) with data and MC uncertainties summed in quadrature. The resulting upper limit of the IC contribution for the Sherpa calculation at NLO w u. l. = 1.93 % is presented in Fig. 1 (bottom). This value of w u. l. corresponds to the χ 2 at minimum plus one, see the solid line in Fig. 3. The χ 2 at forward rapidity is constructed from the data points in bins 3-11. In the χ 2 calculation we include the sum in quadrature of uncertainties from the data (statistical and systematical) and theory (mainly scale uncertainties).
The w cc extraction method was also repeated in the combined QCD approach. The CTEQ66 and CTEQ66c PDFs, which include the IC fraction in the proton, were used to calculate the E γ T -spectrum in the forward rapidity region. The E γ T spectra and χ 2 as a function of w, obtained within these approaches, are presented in Figs. 2, 3 (dashed red line) respectively. The upper limit of the IC contribution obtained within the combined QCD approach is about w u. l. = 2.91 %. As it was shown in [31], the combined QCD approach does not include parton showers and hadronization, a contribution which is sizable at E γ T > 100 GeV, where the IC signal could be visible. Therefore, the results obtained from Sherpa (NLO) approach, which include these effects, are more realistic.
The w-dependence of the χ 2 -functions obtained with both the Sherpa and combined QCD approach in the forward region are presented in Fig. 3. By definition, the minimum of the χ 2 -function is reached at a central value w c which corresponds to the best description of the ATLAS data. The application of Sherpa results in w c = 0.00 % while the combined QCD gives us w c = 1.00 %. Figure 3 shows a rather weak χ 2 -sensitivity to the wvalue. This is due to the large experimental and theoretical QCD scale uncertainties, especially at E γ T > 100 GeV (Figs. 1, 2). Since it is not possible to extract the w c -value with a required accuracy (3-5 σ ) we present the relevant upper limit at the 68 % confidence level (CL). Fwd.
(b) w u. l. = 2.91 % Fig. 2 The spectrum of prompt photons as a function of its transverse energy E γ T calculated with the combined QCD analysis, compared to ATLAS data [32]. a Top: the spectrum in the central rapidity region |η γ | < 1.37 and forward 1.56 ≤ |η γ | < 2.37 region without the IC contribution; a middle: the ratio of the MC calculation to the data for the central rapidity region (w c = 0%); a bottom: the ratio of the MC calculation to the data for the forward rapidity regions (w c = 0%). b The same spectra, as in (a), but with upper limit of IC contribution w u. l. = 2.91 % corresponding to the best fit of the data For an estimation of the scale uncertainty we have used the conventional procedure used in a literature, which is to vary both the renormalization μ R and factorization μ F scales within 0.5E γ T and 2E γ T . In fact, there are several methods to check the sensitivity of observables to the scale uncertainty, see for example [48] and references therein. The renormalization scale uncertainty of the E γ T -spectra poses a serious theoretical problem for obtaining a more precise estimate of the IC probability from the LHC data. Contrary to the NNPDF3.1 estimation of w cc = 0.34 ± 0.14% [26] we are able to estimate only the upper limit of w cc extracted from the ATLAS data [32], as it was shown above. The precision of our analysis is limited by the experimental systematic uncertainties -mainly, by the c-tagging uncertainty which is predominantly connected with the light jet scaling factors [32]. It is also limited by theoretical QCD scale uncertainties. The PDF uncertainties are included in the predictions using the Sherpa (NLO). In contrast to these uncertainties, the statistical uncertainty does not play a large role. In Fig. 4 we have shown the fraction of the uncertainty of the IC upper limit w u. l. at 68 % CL introduced from the uncertainty of each contributing component. The allowed upper limit is presented four times, every time the component of uncertainty in question is reduced from its actual value (100 %). This assumes that the central values of the experiment does not change. For example, according to Fig. 4, if the experimental systematic uncertainty (dashed blue line) is decreased by a factor of 2, then w u. l. is also reduced twice and it becomes about 1 %. One can see also that the IC upper limit is effectively insensitive to a change of the experimental statistical uncertainty. In order to obtain more reliable information on the IC probability in the proton from future LHC data at √ s = 13 TeV, a more realistic estimate of the theoretical scale uncertainties and reduction of the systematic uncertainties is needed.
This problem can be, in fact, eliminated by employing the "principle of maximum conformality" (PMC) [49,50] which sets the renormalization scale by shifting the β terms in the perturbative series of the QCD running coupling. The PMC predictions are independent of the choice of renormalization scheme -a key requirement of the renormalization group. Its utilization will be the next step of our study.
Let us note that the sensitivity of our results to different PDFs is less than the experimental and theoretical scale uncertainties. As shown in Figs. 5, 6 , the difference between E γ T spectra using different PDFs is less than 10% in the central region and less than 5% in the forward rapidity region. In  contrast, the experimental uncertainty is about 50-60%, as it is shown above. Therefore, one can neglect the sensitivity of the E γ T spectrum to PDF uncertainties, especially in the forward rapidity region.

Summary
A first estimate of the intrinsic charm probability in the proton has been carried out utilizing recent ATLAS data of prompt photon production accompanied by a c-jet at √ s = 8 TeV [32]. We estimate the upper limit of the IC prob-   Fig. 6 Top: The E γ T spectrum of photon in the forward rapidity region and with zero IC contribution using different PDFs, namely, CT14nnloIC (black points), CTEQ66 (blue points), CT10 (red points), NNPDF (lilac points). Bottom: the ratio of different PDFs to the CT14nnloIC central set ability in the proton at 1.93% with 68% CL. In order to obtain more precise results on the intrinsic charm contribution one needs additional data and at the same time reduced systematic uncertainties which come primarily from c-jet tagging. In particular, measurements of cross sections of γ + c and γ + b production in pp-collisions at √ s = 13 TeV at high transverse momentum with high statistics [29] will be useful since the ratio of cross sections γ + c/γ + b is expected to be more sensitive to the IC signal [29,30]. The ratio falls to unity slowly with increasing E γ T , in the presence of intrinsic charm. Furthermore, measurements of Z /W + c/b production in pp collision at 13 TeV could also give additional significant information on the intrinsic charm contribution [28][29][30][31]. Our study shows that the most important source of theoretical uncertainty on w cc , from the theory point of view, is the dependence on the renormalization and factorization scales. This can be reduced by the application of the Principle of Maximum Conformality (PMC), which produces schemeindependent results, as well the calculation of the NNLO pQCD contributions. Data at different energies at the LHC which constrain scaling predictions and future improvements in the accuracy of flavor tagging will be important. These advances, together with a larger data sample (more than 100 fb −1 ) at the LHC energy 13 TeVshould provide definitive information on the non-perturbative intrinsic heavy quark contributions and their impact on the fundamental structure of the proton. and less than 5% in the forward rapidity region. The size of the difference is comparable to the size of uncertainty coming directly from the CT14nnlo PDFs and is much smaller than the systematic uncertainty of the data [32].