CT14 Intrinsic Charm Parton Distribution Functions from CTEQ-TEA Global Analysis

We investigate the possibility of a (sizable) nonperturbative contribution to the charm parton distribution function (PDF) in a nucleon, theoretical issues arising in its interpretation, and its potential impact on LHC scattering processes. The"fitted charm"PDF obtained in various QCD analyses contains a process-dependent component that is partly traced to power-suppressed radiative contributions in DIS and is generally different at the LHC. We discuss separation of the universal component of the nonperturbative charm from the rest of the radiative contributions and estimate its magnitude in the CT14 global QCD analysis at the next-to-next-to leading order in the QCD coupling strength, including the latest experimental data from HERA and the Large Hadron Collider. Models for the nonperturbative charm PDF are examined as a function of the charm quark mass and other parameters. The prospects for testing these models in the associated production of a Z boson and a charm jet at the LHC are studied under realistic assumptions, including effects of the final-state parton showering.


I. INTRODUCTION: CTEQ DISTRIBUTIONS WITH INTRINSIC CHARM
The principle of the global analysis is to use QCD theory to analyze a broad range of experimental data, including precision data from HERA, the Tevatron, and the Large Hadron Collider (LHC). In particular, theoretical predictions for short-distance scattering processes allow the measurement, within some approximations, of universal parton distribution functions (PDFs) for the proton. These functions can then be used to predict hadronic cross sections in the QCD and electroweak theories, and in beyond-the-standard-model theories.
With the new high-precision data becoming available from the LHC, the ultimate goal for the global QCD analysis is to be able to make predictions that are accurate to about one percent. This, in turn, requires improvements in theoretical predictions to allow for an accurate extraction of the parton content of the proton in global fits.
A recently published CTEQ-TEA (CT) analysis of QCD data [1] produced the CT14NNLO PDFs, referred to as the CT14 PDFs in this paper. The analysis is based on the next-to-next-to-leading order (NNLO) approximation for perturbative QCD. That is, NNLO expressions are used for the running coupling α S (Q), for the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) evolution equations [2][3][4][5][6], and for those hard matrix elements for which the NNLO approximation is available, such as the deep-inelastic scattering (DIS) neutral-current data from HERA and fixed-target experiments, and the Drell-Yan data from the Tevatron, fixed-target experiments, and the LHC [7][8][9][10][11][12][13][14][15][16][17]. Next-to-leading order (NLO) is used only for inclusive jet data from the Tevatron and the LHC and for deepinelastic scattering (DIS) charged-current data from HERA and fixed-target experiments. The NNLO predictions for these processes [18][19][20] were not available or incomplete at the time of the CT14 study, and we have argued [1,21] that the effect of missing NNLO terms in jet production on the PDFs is small relatively to the experimental uncertainties in the CT14 data sets. Similarly, the NNLO contribution for charged-current DIS, including massive charm scattering contributions, is modest compared to the experimental uncertainties.
In the global analysis, all QCD parameters, such as α s and the quark masses, are correlated with the PDFs. The determination of the PDFs depends not only on the data sample included in the fits, but also on the specific theory assumptions and underlying physics models. As one such choice made in the standard CT PDF sets, the charm quark and antiquark PDFs are taken to be zero below a low energy scale Q c = Q 0 of order of the charm mass. In the CT14 analysis, the charm quark and antiquark PDFs were turned on at the scale Q c = Q 0 = m c = 1.3 GeV, with an initial O(α 2 s ) distribution consistent with NNLO matching [15,22] to the three-flavor result. At higher Q, most of the charm PDF is generated from the DGLAP evolution that proceeds through perturbative splittings of gluons and light-flavor quarks. Hence, the charm PDF from a standard global analysis is called "perturbative", for it was obtained by perturbative relations from light-parton PDFs at scale Q c and perturbatively evolved to the experimental data scale Q.
In addition to the perturbative charm production mechanism, it is believed that "intrinsic charm quarks" may emerge from the nonperturbative structure of the hadronic bound state. The plausibility of the intrinsic charm (IC) component, its dynamical origin, and its actual magnitude have been a subject of a long-standing debate. Indeed, QCD theory rigorously predicts existence of power-suppressed (higher-twist) channels for charm quark production that are independent of the leading-power (twist-2, or perturbative) production of charm quarks. The intrinsic charm (IC) quarks have been associated with the higher |uudcc Fock state of the proton wave function [23][24][25][26][27][28] and predicted by meson-baryon models [29][30][31][32].
On the other hand, Refs. [33,34] concluded that the momentum fraction carried by intrinsic charm quarks is at most 0.5% at the 4σ level, though this conclusion has been challenged in Ref. [35]. This is to be compared to the earlier CT10IC study [36], which concluded that the existing data may tolerate a much larger momentum fraction carried by intrinsic charm quarks. For a valence-like model, it was found to be less than about 2.5%, at the 90% confidence level (C.L.). Recently, several analyses by the NNPDF group [37][38][39][40] established a smaller fitted charm momentum fraction. NNPDF determined a fitted charm momentum fraction equal to (0.26 ± 0.42)% at 68% C.L. just above the charm mass threshold, with the charm quark pole mass taken to be 1.51 GeV [40], and equal to (0.34 ± 0.14)% when the EMC data [41] on SIDIS charm production were included. The current paper revisits the issue in the context of the CT14 analysis [1], also including more recent advances that were made in the follow-up CT14HERA2 study [42]. It updates the previous work [36] on fitting the charm PDFs based on the CT10 NNLO framework [21], as well as the CTEQ6.6 IC study [43] done at NLO. In addition to implementing the combined HERA I+II data on DIS, the new LHC data, and improved parametrizations for light-parton distributions, we shall address some fundamental questions: What dynamics produces the nonperturbative c and c components of the proton? Is there a universal description of this type of charm component that is supported by the QCD factorization theorem, such that the same charm PDF can be used in both lepton-hadron and hadron-hadron scattering processes?
These core questions must be raised to appraise the range of validity of the PDF models with nonperturbative charm in our work and in the other recent studies [33,34,39,40,44]. We address them by starting from the fundamental QCD result, the factorization theorem for DIS cross sections with massive fermions. We start by discussing the definition to the "intrinsic charm", the term that has been used inconsistently in the literature. In the theoretical section, we advance a viewpoint that the "intrinsic charm" can refer to related, but non-equivalent concepts of either the "fitted charm" PDF parametrization, on one hand, or the genuine nonperturbative charm contribution defined by the means of power counting of radiative contributions to DIS. This means that the generic notion of the "intrinsic charm" may cover several kinds of unalike radiative contributions. After we draw this consequential distinction, and assuming that the nonperturbative charm scattering cross section can be approximated by a factorized form, our global analysis examines agreement of various models for the nonperturbative charm with the modern QCD experimental data.
The nonperturbative charm content is normally assumed to be suppressed by powers of (Λ 2 /m 2 c ), where Λ is a nonperturbative QCD scale. But, since this ratio is not very small, it may be relevant in some processes such as precise DIS. The allowed magnitude of the nonperturbative charm is influenced by other theoretical assumptions that a global fit makes, especially by the heavy-quark factorization scheme [22,[45][46][47][48][49], the α s order of the calculation, the assumed charm mass m c , and the parametrization forms for the PDFs of all flavors. We study such effects in turn and find that, among the listed factors, the IC component is strongly correlated with the assumed charm mass.
Dependence on m c in the absence of the nonperturbative charm has been addressed at NNLO in the CT10 NNLO framework [50] and in other references [39,[51][52][53][54][55]. In the context of the CT10 analysis [50], the general dependence on the charm quark mass was studied, and a preferred value of m c (m c ) = 1.15 +0. 18 −0.12 GeV was obtained at 68% C.L., where the error is a sum in quadrature of PDF and theoretical uncertainties. Here, m c (m c ) denotes the running mass of the charm quark, defined in the modified minimal-subtraction (MS) scheme and evaluated at the scale of m c . This value, constrained primarily by a combination of inclusive and charm production measurements in HERA deep-inelastic scattering, translates into the pole mass m pole c = 1.31 +0. 19 −0.13 GeV and 1.54 +0.18 −0.12 GeV when using the conversion formula in Eq. (17) of Ref. [56] at the one-loop and two-loop order, respectively. As the pole mass of 1.3-1.8 GeV borders the nonperturbative region, accuracy of its determination is limited by significant radiative contributions associated with renormalons [57][58][59]. In this light both converted values are compatible with the value of m pole c = 1.3 GeV, which was assumed by CT10 and CT14 and provides the best fit to HERAI+II data at NNLO with the chosen PDF parametric form. We shall use it as our standard charm quark pole mass value in this paper, unless specified otherwise.
To establish robustness of our conclusions, in our fits we varied the selection of data and the analysis setup. Constraints on the IC from both CT14 [1] and CT14HERA2 sets [42] of experimental data were compared. As the CT14HERA2 fit prefers a smaller strangeness PDF than CT14, comparison of the CT14 and CT14HERA2 allowed us to estimate the sensitivity of the IC to the strangeness content. [The sensitivity to the treatment of bottom quarks is expected to be marginal.] Finally, we consider the impact of the possible nonperturbative charm on predictions for the present and future experimental data. The momentum sum rule, one of the key QCD constraints, implies that introduction of a fitted charm PDF modifies the gluon and sea (anti)quark PDFs, particularly, forū andd. Hence, accurate predictions of the c and c parton distributions will be relevant to various important LHC measurements, such as production of W ± , Z 0 , and Higgs boson, or associated production of a charm jet and a Z 0 .
The remainder of this paper is organized as follows. In Sec. II we review the theoretical foundations of the CTEQ global PDF analysis with contributions of massive quarks. In particular, we discuss issues related to the factorization of the charm PDF in the proton, after clarifying the meaning of the PDFs for the leading-power (perturbative) charm, power-suppressed charm, and the fitted charm. Several theoretical models of the intrinsic charm PDF at the Q 0 scale will be presented in Sec. III. The results of our global fits, called the CT14IC PDFs, are discussed in Sec. IV, where the quality of the data description is documented, and a detailed comparison of the CT14IC PDFs with the CT14 PDFs and other PDF sets is provided. The dependence of the CT14IC PDF fits on the charm-quark mass is detailed in Sec. IV C. In Sec. IV F, we discuss the impact of including the EMC data in the global fits for the fitted charm PDFs, as predicted by those theoretical models introduced in Sec. III. We examine the impact of the CT14IC PDFs on the production of the electroweak W ± , Z and Higgs bosons at the LHC in Sec. V, and on a charm jet production associated with a Z boson at the LHC in Sec. VI. Finally, our conclusions are presented in Sec. VII.

BUTIONS
Particle interactions with energies of hundreds of GeV, at modern colliders such as the LHC or the Tevatron, are not directly sensitive to the masses of most Standard Model (SM) fermions. At such high energy, one may safely neglect the mass of any quark in a shortdistance scattering cross section, except for the top quark. Protons, the initial-state nucleons at the LHC, behave as bound states composed of strongly interacting constituents lighter than the top, including light quarks (u, d, s), heavy quarks (c and b), and gluons g. 1 A parton a knocked out of an initial-state proton by a hard collision moves essentially as a massless particle; however, the probability for knocking the parton out, quantified by the parton distribution function f a/p (ξ, µ), or a(ξ, µ) for short, depends on the parton's flavor and, ultimately, the parton's mass.
A charm quark with mass m c ∼ 1.3 − 1.6 GeV is heavier than a proton at rest, with mass 0.938 GeV. If we introduce a parton distribution for the charm, what is the physical origin of this PDF?
The answer is not as clear-cut as for the lighter quarks, whose PDFs are dominated by nonperturbative QCD contributions arising from energies smaller than the proton mass. The light-quark PDFs are essentially nonperturbative; we parametrize each light-quark PDF by a phenomenological function f a/p (x, Q 0 ) at an initial energy scale Q 0 of order 1 GeV and evolve the PDFs to higher energies using the DGLAP equations [2][3][4][5][6]. For the charm and anticharm contributions, on the other hand, the respective PDFs at such low Q 0 are not mandatory. Only some QCD factorization schemes introduce them, with the goal to improve perturbative convergence at scales Q much larger than Q 0 . The perturbative component of the charm PDF dominates in conventional treatments, such as those implemented in the general-purpose QCD analyses by CTEQ-TEA and other groups. However, a nonperturbative component in the charm PDF cannot be excluded either -we will explore it in this paper. What are the theoretical motivation and experimental constraints for the nonperturbative component? Can it be relevant for the LHC applications?
We can systematically approach these questions by reviewing QCD factorization, and the associated factorization theorem, for a perturbative QCD calculation of a radiative contribution with heavy quarks. Let us focus on predictions for neutral-current DIS structure functions F (x, Q) with 3 and 4 active flavors, at a relatively low momentum transfer Q that is comparable to the mass m c of the charm quark. Our considerations can be extended readily to situations with more than four active flavors, and to higher Q values. Moreover, among the experimental processes included in the global QCD analysis, the neutral-current DIS is the most sensitive to charm scattering dynamics [39,50,[52][53][54][55] with the rest of the processes providing weaker constraints. Therefore, it is natural to focus on DIS as the starting point.

A. Exact and approximate factorization formulas
We first write down a phenomenological form for the DIS structure function that is implemented in the CTEQ-TEA PDF analysis: This is a standard convolution formula, consisting of the coefficient function C (N ord ) a (x/ξ, Q/µ, m c /µ; α s (µ)) and the PDFs f (ξ, µ) dependent on the light-cone partonic momentum fraction ξ and factorization scale µ of order Q (set to coincide with the renormalization scale to simplify the notation). The index a denotes the initial-state parton's flavor, running from a = 0, corresponding to the gluon, to the number N f of active quark flavors assumed in the QCD coupling strength α s (µ) and the PDFs f a/p (ξ, µ). Implicitly, summation over quarks and antiquarks is assumed. We reserve the index "h" for a heavyquark flavor, h = c in DIS charm production. 2 The superscripts (N ord ) in both C (N ord ) and f (N ord ) a/p emphasize that their perturbative coefficients are computed up to a fixed order N ord of α s .
Let us highlight several aspects of this formula. First, N f , the number of active flavors, is not measurable, it is a theoretical parameter of the renormalization and factorization schemes chosen for the perturbative calculation. N f should be distinguished from N f s f [60,62], the 2 Beyond the NNLO accuracy considered in this paper, DIS includes contributions with both c and b quarks.
Treatment of such contributions in the ACOT formalism is explained in Refs. [60,61].
number of (anti-)quark species that can be physically produced in the final state in DIS at given collision energy. The optimal value of N f is chosen as a part of the QCD factorization scheme to optimize perturbative convergence. N f s f can be determined from an experimental observable, such as the final-state hadronic mass in the neutral-current DIS process.
Second, the CTEQ-TEA group computes the perturbative coefficients of C in the S-ACOT-χ scheme [46,[63][64][65], a general-purpose factorization scheme for lepton-hadron and hadron-hadron scattering processes. For neutral-current DIS, C (N ord ) a were derived in this scheme up to O(α 2 s ), or NNLO [60]. Figure 1 is reproduced here from Ref. [60] and shows the Feynman diagrams and notations for the perturbative coefficients of the "charm production" structure function F c (x, Q) up to NNLO in the S-ACOT-χ approach. Our discussion will turn to these diagrams for an illustration. The remaining NNLO charm scattering contributions in NC DIS, arising in the light-quark structure function F l (x, Q) and not as important numerically, can also be found in Ref. [60].
Third, in a general-purpose analysis such as CT14 NNLO, we start with non-zero PDF parametrizations for the gluon and 3 light (anti-)quark flavors at the initial scale slightly below the charm mass, Q 0 = m c − ǫ. The input charm mass can be either the MS mass m c (m c ), or the pole mass m pole c : the two are related by NNLO perturbative relations [56,66], both are implemented in CT14 PDFs. 3 As f (N ord ) a/p (ξ, Q) are evolved upward from the initial scale Q 0 , they are converted from N f = 3 to 4, and from 4 to 5, at the corresponding switching points Q i . The perturbative coefficients of C (N ord ) a are converted concurrently to preserve the factorization scheme invariance at each order of α s . The CT14 analysis switches from N f to N f +1 exactly at the heavy quark mass; so for the charm quark the switching takes place at the energy scale Q c = m c .
In this conventional setup, we assume a zero charm PDF, f c/p (ξ, Q 0 ) = 0, for N f = 3 at the initial scale Q 0 slightly below Q c = m c , and obtain a small non-zero f c/p (ξ, Q c ) for N f = 4 at scale Q c via perturbative matching. Of course, Q c is arbitrary, we could equally choose a Q c value below Q 0 and then expect a non-zero charm PDF also at Q 0 . This alternative suggests the possibility of including a non-zero initial charm PDF parametrization, or the "fitted charm" parametrization, at the initial scale Q 0 that would now correspond to N f = 4. However, if the charm quarks are produced exclusively from perturbative gluon splittings, the dependence on the fitted f c/p (ξ, Q 0 ) cancels up to the higher α s order in the cross section, not the PDF alone. It only makes a difference, compared to the higher-order uncertainty, if another mechanism adds up to perturbative charm-quark production.
To demonstrate this, compare the above approximate fixed-order formula (1), which either includes the fitted charm PDF, or not, to the all-order expression for F (x, Q) with massive quarks that follows from the QCD factorization theorem [63,68]: Eq.
(2) underlies all modern computations for the inclusive DIS observables, in the factorization schemes with fixed or varied N f values. The convolution of C a with f a/p (ξ, µ) in Eq. (2) includes all "leading-power" radiative contributions that do not vanish when the physical scales √ s, Q, m c are much larger than the nonperturbative hadronic scale Λ of order less than 1 GeV. In Eq. (1), as implemented in the fits, this leading-power C a ⊗ f a/p is approximated just up to order N ord .
A (2) h,g c h,h  This means that, in the all-order factorization theorem (2), [C a ⊗ f a/p ](x, Q), the first term on the right-hand side, captures all contributions associated with the leading-power, perturbative, charm production. On the other hand, when a non-zero initial condition for f (N ord ) c/p (ξ, Q 0 ) is introduced in the fitted formula (1), it plays the role of a placeholder for several kinds of missing contributions that appear in the full factorization formula (2), but not in the approximate formula (1). For example, it substitutes in part for the leadingpower perturbative contributions beyond the order N ord . The O(α 2 s ), or NNLO, radiative : Dominant squared leading-power amplitudes in DIS charm production in the Q 2 ≫ m 2 c ≫ Λ 2 limit. Here F is the DIS structure function, C 0 , K 0 and D are two-particle irreducible (2PI) subgraphs, T 0 and T (4) 0 are the twist-2 and 4 target hadron subgraphs, and K (2|4) 0 is the heavy-quark "mixed-twist" 2PI subgraph.
contribution to neutral-current DIS heavy-quark production is large numerically. If a global fit is done at NLO, as in Refs. [33,39,43], it prefers an augmented fitted charm f (N LO) c/p (ξ, Q 0 ) of a certain shape in part to compensate for the missing NNLO DIS Wilson coefficients.
The fitted charm may also absorb part of the last, power-suppressed, term on the righthand side of Eq. (2). The "power counting" analysis of Feynman integrals shows that the ordinary power-suppressed contribution in unpolarized inclusive DIS is proportional to (Λ/Q) n with integer n ≥ 2 ("twist-4", see, e.g., [69,70]). In the DIS scattering of charm quarks, the lowest power-suppressed contribution also includes terms of order Λ 2 /m 2 c [23,24,63]. The latter term clearly does not vanish with increasing Q and, furthermore, at very high Q it is enhanced logarithmically and behaves as (Λ 2 /m 2 c ) ln d (Q 2 /m 2 c ) with d ≥ 0 due to contributions from collinear scattering. The power-suppressed charm contribution, once introduced at low scale Q ∼ m c , will survive to the much higher scales relevant to the LHC.

B. Charm contributions in 3-flavor and 4-flavor schemes
While the complete analysis of the twist-4 contribution is far too extensive, we present a heuristic explanation of its logarithmic growth by following an analogy with the leadingpower, or twist-2, terms [60]. It is useful to compare the relevant Feynman graphs in the N f = 3 factorization scheme, the most appropriate scheme to use in the threshold kinematical region, where Q is comparable to m c , and in the N f = 4 scheme, which is most appropriate at Q 2 ≫ m 2 c , where the charm density has the most physical interpretation. First, recall that in the N f = 3 scheme all subgraphs containing heavy-quark propagators are assigned to the Wilson coefficients C a and not to the PDFs f a/p . Among the leading-power hard-scattering amplitudes in Fig. 1, the only contributions arising in the 3-flavor scheme are those attached to the external gluons and light quarks, denoted by F h,a . The explanation for this is that the N f = 3 scheme applies zero-momentum subtraction to UV singularities with heavy-quark propagators and strongly suppresses highly off-shell charm quark propagators as a consequence of their manifest decoupling. Therefore, the non-negligible Feynman integrals in this scheme contain the charm propagators only in the hard-scattering subgraphs, where the virtualities of all particle momenta are comparable to Q 2 and m 2 c . The nonperturbative subgraphs with virtualities much less than Q 2 contain only light-parton propagators, as those are renormalized in the MS scheme.
A twist-4 N f = 3 hard-scattering matrix element for F h,a can be thought of as a twist-2 N f = 3 hard-scattering matrix element connected to the parent hadron by an additional lightparton propagator at any point in the hard subgraph. Both twist-2 and twist-4 terms with charm take the factorized form illustrated in Fig. 2, while Fig. 3 shows representative twist-4 squared matrix elements obtained after attaching the second initial-state gluon to some of the twist-2 matrix elements in Fig. 1. In the hadronic cross section, every twist-4 hard scattering cross section shown in Fig. 3 is multiplied by a twist-4 (double-parton) nonperturbative function, such as f gg/p (ξ 1 , ξ 2 , µ). Insertion of two QCD vertices suppresses the twist-4 cross section by a power of α s compared to the counterpart twist-2 cross section, while the insertion of two propagators and multiplication by a twist-4 function further suppresses it by a power of Λ 2 /p 2 with p 2 of order Q 2 ∼ m 2 c . At twist-4, we encounter several new nonperturbative functions that are not constrained by the data and obey their own evolution equations at the scale Q [71,72]. The complete analysis of twist-4 is lengthy -we will refer to the vast literature on the subject, including Refs. [69,70,[73][74][75][76][77][78][79][80][81][82][83][84][85].
We further note that, in the limit Q 2 ≫ m 2 c ≫ Λ 2 , the twist-4 charm scattering cross sections contain ladder subgraphs of essentially twist-2 topology. They can be seen in Fig. 2, illustrating a decomposition of the structure function F containing the ladder contributions. D denotes a two-particle irreducible (2PI) part (in the vertical channel) of the structure function F . The first graph on the right-hand side is a generic twist-2 ladder contribution recognized from the calculation of NLO splitting functions in the massless case by Curci, Furmanski, and Petronzio [86]. It is composed of 2PI subgraphs C 0 , K 0 , and T 0 (without an upper index), where C 0 and T 0 are coupled to the virtual photon and target hadron, respectively. The decomposition in terms of D and C * K 0 ... * T 0 for twist-2 also appears in the Collins' proof of QCD factorization for DIS with massive quarks [63].
The ladder graphs are different in the N f = 3 and N f = 4 schemes. The N f = 4 scheme introduces additional terms with heavy quarks that approximate the leading contribution in the Q 2 ≫ m 2 c limit. In Fig. 1, these ladders correspond to the contributions proportional to the "flavor-excitation" Wilson coefficient functions c (k) h,h . Such terms are absent in the N f = 3 scheme, and their purpose is to resum collinear logs ln d (Q 2 /m 2 c ) from higher orders with the help of DGLAP equations. In this case both the light-and heavy-parton subgraphs are renormalized in the MS scheme. Importantly, apart from a finite renormalization of α s , the perturbative expansions of the structure functions in the N f = 3 and N f = 4 schemes are equal up to the first unknown order in α s -the condition that we expect to hold both for the twist-2 and twist-4 heavy-quark contributions.
Next to the twist-2 term in Fig. 2 we show a ladder attached to a twist-4 target subgraph T is connected to T 0 by four propagators, at Q 2 ≫ m 2 c ≫ Λ 2 it scales as Λ 2 /p 2 . Since it includes loop integrals with massive quark propagators 1/(k / − m c ), the momentum scale p can be either Q or m c ; but the Λ 2 /m 2 c term is less suppressed than Λ 2 /Q 2 . [It is crucial that two large QCD scales, m c and Q, are present, in contrast to the massless-quark case.] On the other hand, apart from the replacement of T 0 · K 0 by T , the second ladder has the structure of the first one.

Factorization for twist-2 contributions
We assume that the Feynman diagrams in Fig. 2 are unrenormalized and indicate this by a subscript "0". Ref. [63] shows how to recast the full sum of twist-2 diagrams into a factorized convolution by recursively applying a factorization operator Z and renormalizing the UV singularities. Z is a projection operator that is inserted recursively between the rungs of the ladder diagram, e.g., at the location indicated by the circle markers. The action of the Z operator is to replace the exact ladder graph by a simpler, factorized expression which provides a good approximation to the full graph in the Q 2 ≫ m 2 c limit, and which is valid up to a power-suppressed remainder r. In particular, Z replaces the off-shell intermediate parton propagator at the insertion point by an on-shell external state with zero transverse momentum in the Breit frame. By considering recursive insertions of the Z operators to all orders, one demonstrates factorization for F (x, Q) in either the N f = 3 scheme or the N f = 4 scheme of the Aivasis-Collins-Olness-Tung (ACOT) class [46]. By its construction, the remainder r is of order highest virtuality in T 0 lowest virtuality in C 0 with p = Q or m c . While the Z operation in the MS scheme is uniquely defined for intermediate light states, for a heavy quark, it encounters an additional ambiguity. The projection operator acting upon an intermediate heavy quark, denoted by Z h , may include additional powers of (m 2 c /Q 2 ) that vary among the conventions [60,63]. The ambiguity in Z h gives rise to several versions of the ACOT-like schemes, all equivalent up to a higher order in α s . The form of Z h may be even made dependent on the type and α s order of the scattering contribution: some choices for Z h , such as the one made in the SACOT-χ scheme [60,64,65], simplify perturbative coefficients and enable fast perturbative convergence.
In a practical calculation of a twist-2 cross section illustrated by Fig. 1, the Z operation defines the prescription for constructing the perturbative coefficients C where A (k) a,b (k = 0, 1, 2, . . . ) are perturbative coefficients [15] of operator matrix elements for finding a parton a in a parton b. A perturbative coefficient C i,a of the Wilson coefficient function at a s order k can be found by comparing the perturbative coefficients on the left and right sides of The comparison does not specify the form of the perturbative coefficients c (k) i,h with an initialstate heavy quark; those are specified by Z h at each a s order k and re-used in exactly the same form in all occurrences of c  i,h and not the partonic PDF coefficients A (k) a,b that remain defined in the MS scheme. With such self-consistent definition, the dependence on Z h cancels up to the first unknown order in a s , as it was verified numerically up to O(α 2 s ) in Ref. [60]. 4 Up to NNLO, we use a simplified decomposition of the neutral-current DIS structure function over the quark flavors probed by the virtual photons: where e i is the (anti)quark's electric charge [60]. The SU (N f ) decomposition of the ACOT structure functions for higher orders was derived in Ref. [61].

Factorization of twist-4 contributions: a sketch
Going back to Fig. 2, we recall that, while in the twist-2 factorization formula (3) the K (2|4) 0 · T (4) 0 subgraph is counted as a part of the remainder r ∼ Λ 2 /m 2 c , diagrammatically, it is attached to the upper ladder subgraphs in exactly the same way as the twist-2 K 0 · T 0 . We can treat the sum as a modified target contribution of twist-2, which now includes some power-suppressed correction. The derivation of the factorization for Q 2 ≫ m 2 c can be repeated for N f = 4 as in the previous subsection. The factorized cross section reproduces the structure function up to the terms of order Λ 2 /Q 2 or (Λ 4 /m 4 c ). At the level of individual contributions, the K (2|4) 0 · T (4) 0 target subgraph introduces a nonzero term in the charm PDF at the switching scale from 3 to 4 flavors. We can continue to use the DGLAP equations and the same coefficient functions as in the pure twist-2 case, and the latter are again dependent on the definition of operator Z h (the heavy-quark mass scheme). In particular, if the flavor-excitation coefficient function c The net change does not exceed the total error Λ 2 /Q 2 of the factorized approximation.
This implies that the twist-4 component of the charm PDF is compatible with any available version of the ACOT scheme, the differences between the structure functions in these schemes are of order Λ 2 /Q 2 for twist-4 and even weaker for higher twists. Furthermore, by the structure of the ACOT schemes, the scheme differences cancel order-by-order in α s . Therefore, the claim in Refs. [37][38][39] that the nonperturbative charm is only consistent with the "full" version of the ACOT scheme or its analog schemes, such as the fully massive FONLL scheme, is not correct. In our analysis, it suffices to use the S-ACOT-χ scheme, with or without the power-suppressed component. Since open charm is produced in cc pairs in neutral-current DIS, and not as lone c (anti)quarks, the χ rescaling in the S-ACOT-χ scheme [65], requiring production in pairs only, approximates energy-momentum conservation better than its full ACOT counterpart that also tolerates production of single c orc quarks.
Let us illustrate the calculation of the simplest twist-4 charm contributions on an example of select twist-4 squared amplitudes in Fig. 3. Again, we follow a close analogy to the twist-2 S-ACOT-χ calculation in Sec. II B 1, see also [46] and [60].
The first line in Fig. 3 shows the lowest-order twist-4 contributions of order O(α 2 s ), the remaining lines show some radiative contributions of order O(α 3 s ). As before, a superscript in the parentheses indicates the order k of the perturbative coefficient.
In either the 3-or 4-flavor scheme, we start by computing "flavor production" structure functions F (k) h,ab , such as F (2) h,gg or F  coefficients order-by-order as in Eqs. (5) and (6).
For instance, at order α 2 s , the double-convolution integral F h,gg ⊗ ⊗f gg/p over the gluonpair light-cone momentum fractions ξ 1 and ξ 2 scales as α 2 h,h ⊗ f c/p , added across all Q ≥ Q c in order to obtain a smooth prediction for F (x, Q). 5 The h,h ⊗ f c/p may be of the same order as c h,h ⊗ f c/p at relatively low Q. The remainder is given by C h,gg is found from the comparison of the O(α 2 s ) coefficients in Eq. (6): The Feynman diagram for the "subtraction term" c h,gg is shown in the first row of Fig. 3. It is obtained by inserting Z h into the Feynman graph for F (2) h,gg in order to constrain the momentum of the cut charm propagator to be collinear to that of the target hadron, and to replace the part of the graph for F (2) h,gg above the insertion by a simpler subgraph given by c h,h . Clearly, the remainder is process-dependent.
The next-order contribution F h,gg with an added gluon line develops a logarithmic enhancement at m 2 c ≪ Q 2 , which is resummed as a part of c h,gg stands for the infrared-safe part (with respect to light partons) of F h,gg in the MS scheme [60].
The rest of the coefficient functions can be computed along the same lines.

A. Overview
To recap the previous sections, a non-zero initial condition at Q c for the "intrinsic charm PDF", interpreted in the sense of the "fitted charm", may be used to test for the powersuppressed charm scattering contribution of order O(α 2 s ), of the kind shown in Fig. 3. To be sensitive to these contributions, the twist-2 cross sections must be evaluated at least to NNLO to reduce contamination by the higher-order twist-2 terms. The complete set of power-suppressed massive contributions can be organized according to the method of the ACOT scheme. It is comprised of numerous matrix elements F h,b for double-parton initial-states b, as well as of twist-4 nonperturbative functions such as f gg/p (ξ 1 , ξ 2 , µ).
Various model estimates suggest a power-suppressed charm cross section of a modest size: of order of a fraction of the α 2 s component in DIS charm production, carrying less than about a percent of the proton's momentum. To estimate sensitivity of the QCD data before resorting to the full twist-4 calculation, we utilize an update of the phenomenological method of the CTEQ6.6 IC NLO and CT10 IC NNLO analyses [36,43]. In contrast to the previous analyses, we examine a more extensive list of nonperturbative models, fit the most complete set of DIS data from HERA as well as the data from the LHC and (optionally) the EMC, and utilize a PDF parametrization that results in a more physical behavior.
Four models for the charm-quark PDF c(x, Q 0 ) ≡ c(x) at the initial scale Q 0 will be considered. [Q 0 is set to be less than Q c = m pole c in all cases.] Besides the conventional CT14 model that sets c(x) = 0, the other three models allow for c(x) of an arbitrary magnitude. In all models, the charm PDF is convoluted with the S-ACOT-χ coefficient functions c It remains constant below the switching scale Q c and is combined with the perturbative charm component at Q c and evolved to Q > Q c by the 4-and 5-flavor DGLAP equations.
Neither the present fit, nor the contemporary fits by the other groups include the twist-4 remainders of DIS cross sections discussed in Sec. II B 2: C h,gg ⊗ f gg/p , etc. The remainders are process-dependent and comparable to the c h,h ⊗ f c/p convolutions at energies close to m c . Without including these process-dependent terms explicitly, the fitted charm PDF found in a fit to DIS is not a truly universal nonperturbative function; it absorbs the above process-dependent remainders. Furthermore, in DIS at very low Q or W , separation of the Λ 2 /Q 2 and Λ 2 /m 2 c terms presents an additional challenge. The experimental data in the CT14(HERA2) fits is selected with the cuts Q 2 > 4 GeV 2 , W 2 > 12.5 GeV 2 so as to minimize sensitivity to the Λ 2 /Q 2 terms. This is usually sufficient to minimize the Λ 2 /Q 2 contributions below the PDF uncertainty from other sources. We examine the possibility of the impact of the Λ 2 /Q 2 terms on the best-fit c(x, Q 0 ) in Sec. IV E.

B. Valence-like and sea-like parametrizations
Given that several mechanisms may give rise to the fitted charm, we will parametrize it by two generic shapes, a "valence-like" and a "'sea-like" shape. The two shapes arise in a variety of dynamical models.
A valence-like shape has a local maximum at x above 0.1 and satisfies The distributions for valence u and d quarks fall into this broad category, as well as the "intrinsic" seaquark distributions that can be naturally generated in several ways [25], e.g., for all flavors, nonperturbatively from a |uudQQ Fock state in light-cone [23,24,[26][27][28] and meson-baryon models [29][30][31][32]; forū andd, from connected diagrams in lattice QCD [87]. 6 A sea-like component is usually monotonic in x and satisfies x → 1, with a 1 slightly above 1, and a 2 5. This behavior is typical for the leading-power, or "extrinsic" production. For example, an (anti)quark PDF with this behavior originates from g → qq splittings in perturbative QCD, or from disconnected diagrams in lattice QCD (see Ref. [87] for details). Even a missing nextto-next-to-next-to-leading order (N3LO) leading-power correction may produce a sea-like contribution at x ≪ 0.1, where the valence-like components are suppressed. One may wonder why the charm quark PDF cannot be fitted to a more general parametrization, in the same manner as the light-quark PDFs. We find that the primary problem is that there are not enough precision data available to provide meaningful constraints on the power-suppressed IC content in the {x, Q} regions where it can be important (see the discussion of the EMC charm data in Sec. IV F). There is also a danger that the charm quark distribution, being relatively unconstrained, may behave unphysically, for example, when the fit allows a valence-like c(x, Q 0 ) to be almost the same in size asū(x, Q 0 ) or d(x, Q 0 ) at Q 0 ∼ m c and x → 1, where the experimental constraints are weak. We must also demand conceivable cross sections to be non-negative, even though the PDFs themselves can generally have a negative sign. Adopting a too flexible fitted charm PDF parametrization may mask unrelated higher-order radiative contributions to the data, hence lead to misinterpreted fits. Thus, we restrict the freedom of the charm quark somewhat by constraining it to be non-negative and have either a valence-like or sea-like form, with only one free multiplicative parameter. The positivity of the BHPS form enables positive charm-scattering cross sections at large x, while a negative-valued SEA form is not statistically distinguishable in the fit from a positive SEA form at a larger m c value. [The dependence of SEA fits on m c is reviewed in the next Section.] We have verified that a mixed charm parametrization that interpolates between the valence-like and sea-like parametrizations only slightly increases the range of the allowed charm momentum fraction, without impacting the main outcomes.

C. The charm distribution models in detail
We will now review these four models, whose xc(x, Q) distributions at Q = 1.3 GeV and Q = 2 GeV are depicted in Fig. 4 for later reference. These models are implemented in five fits, BHPS1,2,3 and SEA1,2, summarized in the next section.
i) Perturbative charm. The first model is the one used in the standard CT14 (and CT14HERA2) PDF fits, in which a non-zero charm PDF is produced entirely perturbatively by NNLO switching from the 3-flavor to the 4-flavor scheme at the scale Q c . The size of the preferred charm distribution at a given Q significantly depends on the values of the physical charm quark mass m c and QCD coupling strength α s (m Z ). On the other hand, its dependence on the auxiliary theoretical scales of order m c , including the switching scale Q c and the scale in the rescaling variable χ, cancels up to N3LO and thus is relatively weak; see a practical illustration in Fig. 1 of [50]. The net momentum fraction of the proton carried by charm quark starts off close to zero at Q ≈ Q c and effectively saturates at high Q values at a level of approximately 2.5%, see Fig. 7.
ii ) The approximate Brodsky-Hoyer-Peterson-Sakai (BHPS) model [23,24] parametrizes the charm PDF at Q 0 by a "valence-like" nonperturbative function This function is obtained from a light-cone momentum distribution by taking the charm mass to be much heavier than the masses of the proton and light quarks: Here and in the following, A is the normalization factor that is to be determined from the fit. This parametrization choice is employed in two global fits named BHPS1 and BHPS2, corresponding to two values of A in Eq. 10. The parametrizations forū(x, Q 0 ) andd(x, Q 0 ) in this case are taken to be the same as in the CT14/CT14HERA2 fits, i.e., they do not have a "valence-like" component and monotonically decrease at x → 1. The parametrizations of this kind tend to have enhancedc/ū andc/d ratios at x → 1, see Fig. 9.
iii) The exact solution of the BHPS model is realized in the BHPS3 fit. Instead of approximating the probability integral as in model ii ), the c(x) is obtained by solving the BHPS model for the |uudcc Fock state numerically and keeping the exact dependence on M p , m u , and m d . This fit also includes small BHPS contributions to theū andd antiquarks generated from the |uuduū and |uuddd Fock states according to the same method. In the BHPS model, the quark distributions are determined by starting from a |uudqq proton Fock state, where the probability differential for a quark i to carry a momentum fraction x i is given by The standard BHPS result, used in ii ), is given by letting q = c and taking the limit m c ≫ M p , m u , m d to produce Eq. (10). However, Ref. [26] has shown that the solution that keeps the masses finite, including those of the light quarks, modifies the shape of c(x), slightly shifting the peak to smaller x. A similar conclusion was reached in Ref. [27], where a kinematic condition on the intrinsic charm was determined analytically by neglecting the masses of the three light valence quarks and retaining the ratio M 2 p /m 2 c . The change in the BHPS charm quark PDF from including the full mass dependence, although visible, is small compared to the uncertainties in the global analysis. However, by using this generalized BHPS model (BHPS3) in the context of the CT14HERA2 fit, and also including the BHPSū andd components, we obtain physically consistent ratios of the charm-quark and light-antiquark PDFs at large x, cf. Fig. 9. We do not, however, include the BHPS contribution to the s quark PDF, because it is overwhelmed by the very large strange PDF uncertainty. The presence of a BHPS component for the strange quark does not affect our conclusions about the nonperturbative charm, so we leave this topic for a separate CTEQ study of the strange content of the proton.
iv ) In the SEA model, the charm PDF is parametrized by a "sea-like" nonperturbative function that is proportional to the light quark distributions: This model is assumed with the SEA1 and SEA2 PDF sets from the two global fits distinguished by the value of normalization A in Eq. 12. Finally, the normalization coefficient A in models ii )-iv ) can be derived from the charm momentum fraction (first moment) at scale Q: By its definition, x IC is evaluated at the initial scale Q 0 . It is to be distinguished from the full charm momentum fraction x c+c (Q) at Q > Q c , which rapidly increases with Q because of the admixture of the twist-2 charm component.

A. Settings of the fits
The BHPS1, BHPS2, SEA1, and SEA2 parametrizations are obtained by following the setup of the CT14 analysis [1]. BHPS3 is obtained with the CT14HERA2 setup [42]. The CT14HERA2 NNLO fit is very similar to the CT14 fit except that the HERA Run I and II combined cross sections were used in place of the Run I cross sections. One of the poorly fit NMC data sets [89] was dropped in CT14HERA2, and the low-x behavior of the strange (anti)quarks was no longer tied to that of theū andd antiquarks. This extra flexibility in s(x, Q 0 ) of CT14HERA2 resulted in a reduction of s(x, Q 0 ) over the entire x range relatively to CT14. This feature has potential implications for the models of c(x) with a sea-like behavior. In some exploratory fits, we include the EMC data [41] on semiinclusive DIS charm production, while in the other fits we examine sensitivity on the input pole charm mass. 7 The PDFs for light partons are parametrized at an initial scale slightly below Q 0 = m pole c = 1.3 GeV, with the exception of the study of the m pole c dependence, in which it was more convenient to start at a lower initial scale Q 0 = 1.0 GeV. For all models, the QCD coupling constant is set to α s (M Z ) = 0.118, compatible with the world average value [67] α s (M Z ) = 0.1184 ± 0.0007, as in the standard CT PDF fits. The PDFs are evolved at NNLO with the HOPPET code [90]. NLO ApplGrid [91] and FastNLO [92] interpolation interfaces, combined with NNLO/NLO factor look-up tables, were utilized for fast estimation of some NNLO cross sections.

B. Dependence on the charm momentum fraction
In the models in Sec. III C, the magnitude of c(x) is controlled by its normalization A, correlated uniquely with the net momentum fraction x IC of c(x, Q 0 ) +c(x, Q 0 ) defined in Eq. (13). The choice of the x IC affects theoretical predictions in a number of ways, either directly by modifying the charm scattering contributions, or indirectly via the proton sum rule that changes the momentum fractions available to other parton flavors.
To gauge the preference of the global QCD data to a specific x IC , we examine the goodness-of-fit function constructed in the CT14 method from the global χ 2 global and a "tier-2" statistical penalty P [1]. It is convenient to compare each fit with an x IC = 0 to the "null-hypothesis" fit obtained assuming x IC = 0. Thus, we start by computing where χ 2 and χ 2 0 are given for x IC = 0 and x IC = 0, respectively, at 50 values of x IC and default Q 0 = m pole c = 1.3 GeV. We plot the resulting ∆χ 2 behavior in Fig. 5. The CT14 (CT14HERA2) data sets are compared against the approximate (exact) solution of the BHPS model, respectively. The SEA charm parametrizations are constructed as in Eq. (12) in terms of the respective CT14 or CT14HERA2 light-antiquark parametrizations.
For each series of fits, we show curves for two types of estimators: a dashed curve for ∆χ 2 global without the tier-2 penalty P , and a solid one for χ 2 that includes P according to Eq. (14). The χ 2 global function estimates the global quality of fit and is equal to the sum of χ 2 contributions from all experiments and theoretical constraints. A non-negative "Tier-2" penalty P is added to χ 2 global to quantify agreement with each individual experiment [21,36]. Being negligible in good fits, P grows very rapidly when some experiment turns out to be inconsistent with theory. The net effect of P is to quickly increase the full χ 2 if an inconsistency with some experiment occurs, even when χ 2 global remains within the tolerable limits.
We see from Fig. 5 that large amounts of intrinsic charm are disfavored for all models under scrutiny. A mild reduction in χ 2 , however, is observed for the BHPS fits, roughly at x IC = 1%, both in the CT14 and CT14HERA2 frameworks.
The significance of this reduction and the upper limit on x IC depends on the assumed criterion. In CTEQ practice, a set of PDFs with ∆χ 2 smaller (larger) than 100 units is deemed to be accepted (disfavored) at about 90% C.L. Thus, a reduction of χ 2 by less than forty units for the BHPS curves has significance roughly of order one standard deviation. We also obtain the new upper limits on x IC in the CT14 and CT14HERA2 analyses at the 90% C.L.: x IC 0.024 for CT14HERA2 BHPS, x IC 0.016 for CT14 and CT14HERA2 SEA.
In keeping with the previous analysis of Ref. [36], we define specific fits with particular choices of x IC for both examined models. The fits BHPS1 and SEA1 correspond to x IC = 0.6%, while BHPS2 has x IC = 2.1% and SEA2 has x IC = 1.6%. Both the BHPS2 and SEA2 charm parametrizations lie near the edge of disagreement with some experiments in the global analysis data according to the CTEQ-TEA tolerance criterion, cf. Fig. 5. In the CT14HERA2 fit, the BHPS3 point corresponds to x IC = 1%, which represents the best-fit momentum fraction in the CT14HERA2 analysis. We remind the reader that, in addition to fitting more recent experimental data from the LHC and other experiments, the BHPS3 analysis also employs a general numerical solution to the BHPS probability distributions and small valence-like contributions for both theū andd quarks.
The results in Fig. 5 are compatible with the findings of the previous CT10 NNLO IC analysis [36]. In particular, comparing to CT14 in the left frame of Fig. 5 and to Fig. 2 in Ref. [36], 8 we see that the minimum in ∆χ 2 in the right frame of Fig. 5 deepened by approximately 10 units for BHPS3/CT14HERA2 -a minor reduction caused mostly by the change to the CT14HERA2 setup, either for the exactly or approximately solved BHPS model.
Also, for the CT14HERA2 analysis in Fig. 5 (right), we note that ∆χ 2 of the SEA model rises more rapidly with increasing x IC than it does in the comparable CT14 fit. This is due to the greater flexibility in the low-x behavior of the strange-quark distribution in the CT14HERA2 framework discussed previously. More freedom reduces s(x, Q) at low x and thus increasesū(x, Q) andd(x, Q) at the same x. In the CT14 fit with the SEA charm component, the ∆χ 2 minimum is at x IC ≈ 0.004, and it is largely washed out in the CT14HERA2 case. The ∆χ 2 for SEA grows faster for CT14HERA2 compared to CT14: at x IC = 1.6% it is higher by about 40 units in Fig. 5(right) relatively to Fig. 5(left).
The reduction in χ 2 for the NNLO BHPS fits at x IC = 0.01, relatively to the fit with x IC = 0, thus remains a persistent feature of the CT10, CT14, and CT14HERA2 analyses. While the ∆χ 2 reduction is not statistically significant, it raises one's curiosity: is it a sign of a genuine charm component or of the other circumstantial factors identified in Sec. III A? It will be discussed in Sec. IV E that χ 2 is reduced primarily in a few fixed-target experiments (the F 2 measurements from BCDMS and the E605 Drell-Yan data) that are not overtly sensitive to charm production. Conversely, the description of the other experiments that might be expected to be most sensitive to intrinsic charm is not improved.

C. Dependence on the charm-quark mass and energy scale
We have checked that these conclusions are not strongly dependent on the PDF parametrizations of the light partons. However, the SEA parametrization at the initial Q 0 is very sensitive to the assumed charm mass.
Distinct from the auxiliary QCD mass parameters -Q 0 , Q c , and the mass in the χ rescaling variable -the physical charm-quark mass of the QCD Lagrangian enters the DIS hard matrix elements through the "flavor-creation" coefficient functions, such as the ones for the photon-gluon fusion. The NNLO fit to DIS is mostly sensitive to the primordial QCD mass parameter m c , not to the auxiliary parameters of order m c [50]. The m pole c dependence 8 In Ref. [36], χ 2 global , P , and χ 2 are denoted by χ 2 F , T 2 , and χ 2 F + T 2 .  Fig. 6. The general setup of these χ 2 scans follows the fits to the CT14 data. To access the masses below 1.3 GeV, we reduced the initial scale Q 0 to 1 GeV and examined alternative forms for the gluon PDF parametrization, because DIS charm production is sensitive to the gluon PDF g(x, Q). Dependence of ∆χ 2 for CT14 NNLO on m pole c for two representative gluon parametrizations at Q 0 = 1 GeV, dubbed "gluon 1" and "gluon 2", is shown in the upper inset of Fig. 6. With the "gluon 1" parametrization, used in the default CT14 fit with Q 0 = 1.3 GeV, g(x, Q 0 ) is constrained to be positive at all x; while for "gluon 2", it is allowed to be negative at the smallest x and Q, provided that the negative gluon does not lead to unphysical predictions. In the latter case, an additional theoretical constraint was enforced to ensure positivity of the longitudinal structure function F L (x, Q) measured by the H1 Collaboration [93]. The more flexible "gluon 2" parametrization results in a marginally better χ 2 with respect to the nominal CT14, or "gluon 1", at a slightly lower m pole c = 1.22 GeV, and with a large uncertainty. This best-fit m pole c value in this range is consistent with the previously observed tendency of the DIS data to prefer smaller MS masses at O(α 2 s ), e.g., m c (m c ) = 1.15 +0.18 −0.12 GeV obtained in the CT10 setup [50].
The two lower insets of Fig. 6  increases, the twist-2 γ * g fusion contribution in the inclusive DIS structure functions is reduced due to phase-space suppression. This suppression is compensated by allowing a larger magnitude of intrinsic c(x, Q), which enhances the γ * c scattering contribution. An opposite effect occurs when m c decreases (i.e., less phase-space suppression for γ * g fusion, a smaller intrinsic charm momentum allowance in γ * c scattering). But theū andd quark PDFs are well constrained by the data, especially from novel cross section measurements for vector boson production in pp and pp in the intermediate/small x region. The net effect is the ∆χ 2 enhancement in the sea-like scenario for m pole c < 1.3 GeV and also for larger x IC fractions. To conclude the discussion of the partonic momentum fractions, Fig. 7 illustrates the first moments x (Q) of the other parton flavors as a function of the factorization scale Q. The momentum fractions are computed separately for quarks, antiquarks, and gluons in the context of the CT14 setup. In the two upper subfigures, the PDF first moments are shown for the BHPS model, while those from the SEA model are shown in the lower two subfigures. The dashed curves represents BHPS1 (SEA1), the dotted ones represent BHPS2 (SEA2).
The lower part of each figure shows x normalized to its CT14 central value. The BHPS2 model curve lies on the edge of the allowed CT14 u and d quark uncertainties, while the SEA2 is on the boundary of theū andd uncertainties. This corroborates the earlier statement that BHPS2 and SEA2 are the extreme choices for the valence-like and sea-like charm distributions, respectively. Next, we will consider the full x dependence of the PDFs provided by our models.

D. Impact of IC on the PDFs
To complement the visualization in Fig. 4 of the x dependence of the BHPS/SEA charm quark PDFs, in Fig. 8 these PDFs are shown normalized to the CT14 charm PDF with no IC contribution. The blue shaded region represents the CT14 uncertainty for c(x, Q) at the 90% C.L.
At low scales (Q = 2 GeV), the charm quark in the SEA models, especially the SEA2 model, appears to be larger, with respect to the CT14 central charm, over a wide range of momentum fraction x. The charm quark distributions from both of these models are clearly outside the CT14 uncertainty bands. Of course, this is not a contradiction, since the CT14 charm PDF is purely radiative, and so it depends on the theoretical assumptions in addition to the constraints from the experimental data. The inclusion of nonperturbative sources of charm relaxes the theoretical assumptions, and so allows a larger charm PDF. The SEA models exhibit minor shape distortions; two bumps are present in both the SEA1 and SEA2 models at x ≈ 10 −3 and x = 0.1.
The charm-quark distributions in the BHPS models at low scales are basically coincident with CT14 below x ≈ 5 × 10 −2 , while a rapid growth is observed at high x, of the largest rate for the BHPS2 model. We note that there is no qualitative difference in the behavior of c(x, Q) between the BHPS3 model and the other BHPS models below x ≈ 5×10 −2 , while the differences at larger x can be ascribed to the exact solution for mass dependence in BHPS3. At a higher scale (Q = 100 GeV), the excesses for all models are suppressed for x 10 −2 due to the effects of DGLAP evolution. The results for the ratio of c(x, Q) IC /c(x, Q) CT14HERA2 are analogous to those shown in Fig. 8 and are omitted.
Additional insights can be gathered by examining the ratios of the charm-quark PDF to other flavors: (c(x, Q) +c(x, Q)) / ū(x, Q) +d(x, Q) , c(x, Q)/u(x, Q), and c(x, Q)/d(x, Q). These ratios are plotted versus x in Fig. 9, for two different values of the Q scale. Also shown for a comparison are the corresponding CT14 PDF uncertainty bands.
For (c +c)/(ū +d), all the BHPS and SEA models reproduce the shape of CT14 at low x, with the ratios in the SEA models shifted upwards. The SEA models retain the shape of CT14 (but with a larger normalization) at higher x as well. All BHPS ratios start to rise quickly in the range 0.1 x 0.2. This rise is essentially unabated at x > 0.2 for the BHPS1 and BHPS2 models, because their respective parametrizations forū andd fall off as (1 − x) d and are more strongly suppressed at x → 1 than the BHPS charm quark PDF. Inclusion of the intrinsicū andd components in the BHPS3 model, together with the numerical estimation of the BHPS integrals for thec,ū, andd intrinsic parametrizations, results in a softer BHPS3 (c +c)/(ū +d) ratio at large x with a bump residing at x ≈ 0.5. The exact amount of suppression at x > 0.5 can be determined, e.g., by a fit to the numerical solutions of the BHPS3 model. In particular, we find that a 6-parameter fit using f (x) ∝ x p 1 (1 − x) p 2 (1 + p 3 x p 4 + p 5 x + p 6 x ln (x)), gives a large-x suppression power p 2 ≈ 8, 9, 10 for intrinsicc,d, andū, respectively.
The c(x, Q)/u(x, Q) ratios in all BHPS models agree with CT14 over the range 10 −5 x 0.1 and exhibit a bump (most prominent for BHPS2) at x ≈ 0.5. The SEA model ratios are notably larger than CT14 in the range 10 −5 x 0.3 and approach CT14 for larger x-values. At higher scale, Q = 100 GeV, all models are closer to CT14 over the range 10 −5 x 0.1 with the exception of SEA2, while the bump in the BHPS models at x ≈ 0.5 are slightly suppressed. The c(x, Q)/d(x, Q) ratio plot shows essentially the same features as the c(x, Q)/u(x, Q) plot, with the difference that the bumps present in the BHPS1, BHPS2 and BHPS3 models, at x ≈ 0.5, are much more pronounced.
An additional charm component (either a sea-like or valence-like one) affects both those LHC predictions that directly involve charm quarks in the initial state, and those that do not. In Fig. 10 we show how the gluon-gluon luminosity is affected by BHPS and SEA models at LHC run I and II energies in the x range sensitive to Higgs production. The parton luminosity is defined as in Ref. [94]. The various models, shown as ratios to CT14NNLO, are well within the 68% C.L. PDF uncertainty. At √ s = 8 TeV the most prominent deviations are for the SEA2 model, which is suppressed at lower M X and is notably larger than CT14 for M X in the TeV range. The BHPS models are almost coincident with CT14 for the invariant mass M X < 200 GeV: BHPS1 and BHPS2 are highly suppressed above M X > 300 GeV, while BHPS3 is suppressed for 0.3 < M X < 3 TeV and enhanced above this energy by approximately 3%. The impact on the Higgs cross section is small, the influence on the high-mass gg PDF luminosities is more pronounced, but still within uncertainties.

E. Agreement with experimental data sets
In this section we focus on the data sets whose goodness-of-fit values are affected by the introduction of the intrinsic charm component. These are selected by computing an effective Gaussian variable, S n , for each experiment n, according to the method introduced in Refs. [21,36,95]. For specifications of S n , we refer the reader to the appendix of Ref. [36]. S n maps the goodness-of-fit χ 2 n for a particular data set, assumed to obey the chi-square probability distribution with N pts data points, onto a variable S n , which obeys a standard normal distribution independently of N pts . More precisely, S n is defined so that the cumulative standard normal distribution evaluated at S n equals the cumulative χ 2 (χ 2 n , N pts ) distribution evaluated at χ 2 n . We adopt an accurate approximation for S n given by The S n distribution over the individual data set characterizes the agreement with the totality of the fitted experiments, regardless of their numbers of data points. Conversely, a naive use of the global χ 2 as the only discriminating variable may give too much weight to the data sets with large numbers of data points, even if the correlations with the fitting parameters are not very significant. The values of S n can easily be interpreted in terms of the probabilities associated with a normal distribution. Fits with S n between -1 and 1 are accepted as reasonable, within the 68% C.L. uncertainties. That is, an increase of S n by 1 has about the same significance (68%) as the increase of χ 2 n /N pts by 2/N pts . Fits with S n > 3 are considered poor, while those with S n < −3 actually fit the data much better than one would expect from the regular statistical analysis: for some reason they have anomalously small residuals.
In Fig. 11, we selectively plot S n for those data sets whose agreement with theory is most affected by the IC in the CT14 fit with m pole c = 1.3 GeV. S n is plotted as a function of x IC for both the BHPS (left) and SEA (right) models.
For the BHPS model, the most visible dependence is found for the fixed target measurements from BCDMS for F p 2 and F d 2 (ID 101, 102) [96,97] and the ATLAS 7 TeV W/Z cross section measurements [98] (ID 268). The E866 Drell-Yan dimuon cross section measurement [99] also shows some variation, however, its S n is always larger than 3 and not shown in Fig. 11(left). These experiments, mostly sensitive to u and d quarks at large x, (slightly) favor a non-zero intrinsic charm component. Although the improvement for the BCDMS S n is relatively mild, the two data sets contain a large number of data points (N pts = 339 for F p 2 and 251 for F d 2 ). The shallow minimum of 20-30 units occurring in χ 2 for the BHPS model in Fig. 5 is attributed primarily to these two experiments; it is not clear whether it originates from the charm component or reflects a small admixture of the N3LO contributions or even some residual 1/Q 2 terms that may be present at relatively low Q and large x.
Continuing with BHPS, the charged-current (CC) DIS measurement [100] F p 2 by CCFR (ID 110) has 0 < S n < 1 for 0 ≤ x IC ≤ 0.02, then S n increases faster for even larger x IC . The combined HERA charm production [101] (ID 147) exhibit 1 < S n < 2 over the whole range of x IC .
The S n dependencies of various experiments for the SEA model are shown in the rightside of Fig. 11. The HERA charm production and BCDMS (F p 2 ) data are very sensitive to x IC in the SEA model. A fast growth for S n is observed for x IC > 0.01, paralleling the increase in χ 2 observed in

F. A global analysis including the EMC charm DIS measurements
The measurement of semi-inclusive dimuon and trimuon production in DIS on an iron target by the the European Muon Collaboration (EMC) [41]   1983, did not follow the stringent criteria on the documentation of systematic uncertainties adopted in more recent studies; therefore, there is a lack of the control on the constraints that these data may impose. This is why the EMC measurements are not included in the CTEQ PDF analyses, whose policy is to include only data with documented systematic errors. Moreover, the EMC analysis has been done at the leading order of QCD, clearly insufficient for accurate conclusions at NNLO. Despite the tensions 9 stated between the EMC measurement and its contemporary experiments in the case of inclusive DIS [89,96,97,105] and semi-inclusive charm DIS production cross sections [41,106] 10 , various studies [107][108][109][110][111] have interpreted the excess seen in a few high-x bins of the EMC F 2c (x, Q) data as evidence for some nonperturbative charm contribution, while yet other studies concluded the opposite [31,33,34]. Our special series of the CT14 IC fits included the EMC F 2c (x, Q) data to investigate the above conclusion. We observe that the EMC F 2c data do not definitively discriminate between the purely perturbative and intrinsic charm models, hence we do not include them in the final CT14 BHPS and SEA fits. However, it is still useful to examine how the EMC data could possibly affect the amount of the intrinsic charm-quark content, especially given their emphasis in a recent NNPDF study [39]. Our findings concerning the fit to the EMC data can be summarized as follows.

χ 2 values for the EMC data set
Either by fitting to the EMC F 2c data or not, we obtain χ 2 /N pts between 2.3 and 3.5 for the EMC data set in various candidate fits. So, for their nominal experimental errors, the EMC data is in general not fit well in either CT14 or CT14HERA2 setup, regardless of the charm model. On the other hand, these χ 2 /N pts values are not dramatically high, it may be argued that allowing for a modest systematic error would improve the agreement to tolerable levels. One way or another, the unknown systematics of this measurement prevents us from concluding for or against the preference of the EMC F 2c data for a particular charm model. To show an example of this, Table I reports the values of χ 2 /N pts for all experiments, HERA inclusive DIS, HERA charm SIDIS, and EMC charm SIDIS in the CT14 (CT14HERA2) NNLO IC candidate fits in the upper (lower) half of the table.
The first two lines in each half present the fits without the nonperturbative charm. When χ 2 for the EMC F 2c data is included with weight 0 (so that the EMC F 2c data has no effect on the PDFs), we obtain χ 2 /N pts ≈ 3.5 -it is quite poor. When the EMC weight is increased to 10 to emphasize its pull, χ 2 /N pts decreases to 2.4, at the cost of a worse χ 2 for the inclusive HERA I+II data and other experiments, and somewhat better χ 2 for charm DIS hadroproduction. Again the quality of the fits is poor, yet it is also compatible with the possibility of moderate unaccounted systematic errors, as those are unknown in the EMC case.
We can also see from Table I that including the BHPS intrinsic charm does not qualitatively change the fit to the EMC data. Without the IC, the χ 2 for all experiments slightly grows if we increase the weight of the EMC data set; with the BHPS intrinsic charm, there seem to be no effect with and without the EMC data, as χ 2 does not change in either case. In the SEA model fit, inclusion of the EMC data results in a larger χ 2 with respect to the fits without the intrinsic charm; description of both HERA inclusive DIS and HERA combined charm SIDIS production deteriorates. To summarize, in all considered intrinsic charm models (BHPS, SEA, and the mixed model that produces a similar outcome), the intrinsic charm has no decisive effect on improving the fit to the EMC data. Figure 12 compares the dependence of ∆χ 2 on x IC in the context of the CT14 and CT14HERA2 global analyses with and without EMC data. It must be noted upfront that, since the EMC F 2c (x, Q) data set are not well described, these ∆χ 2 scans do not establish clear-cut constraints on x IC , contrary to the CT14 IC fits without the EMC data set that were presented earlier.

Constraints from EMC on the IC momentum fraction
The

CT14HERA2, and (b) various gluon PDF parametrizations utilized.
For the BHPS model in Fig. 12 (left), we observe two distinct trends in the fits with and without the EMC data. The spread in the ∆χ 2 band without the EMC data is mostly driven by the differences in the data sets and in the strangeness parametrization between CT14 and CT14HERA2 (the dependence on the gluon parametrization is weak). Meanwhile, after including the EMC data, the spread due to the gluon parametrization dependence is much larger and gives the major contribution to the band. The BHPS model is affected more by the EMC data, the ∆χ 2 band narrows near the minimum when these data are included. The χ 2 minimum with the EMC data moves to a lower value of x IC ≈ 0.006, with substantially the same χ 2 (same depth) at the minimum. The nominal upper limit on x IC moves to about 0.012; its exact location is debatable because of the overall poor quality of the EMC fit, see above.
To contrast with the BHPS case, in the SEA model in Fig. 12 (right), the ∆χ 2 behavior is only mildly impacted by the EMC data. As already discussed in Section IV B and shown in Fig. 5, the ∆χ 2 trend in the SEA model is mostly affected by the differences between the CT14 and CT14HERA2 fits. The EMC data do not change this trend. Both minima are shallow and higher than in the BHPS case.
The Gaussian variables S n quantifying the agreement with the individual data sets are shown for the CT14 fits and for various x IC values in Fig. 13. [The behavior of S n in the CT14HERA2 fit is largely analogous.] In this figure we selected only the experiments that have pronounced dependence on x IC .
Comparing Fig. 13 with Fig. 11 in which the EMC data are not included, one sees that the dependence of S n for the non-EMC experiments on x IC does not qualitatively change upon the inclusion of the EMC. The S n value for the EMC F 2c , indicated as "experiment ID 170", is very high for any x IC . In the BHPS model in Fig. 13 (left), the S n variable for the EMC experiment increases rapidly past x IC of about 0.005, up to very high values at x IC = 0.03. The tier-2 contribution associated with the rapid increase of this S n above 6 produces the rapid rise of the global ∆χ 2 for x IC > 0.01 in Fig. 12. In the SEA model in Fig. 13 (right), we observe S n > 4 for the EMC regardless of x IC . To recap, the EMC data has a weak impact on fitting the rest of the CT14/CT14HERA2 data. Increasing the weight of the EMC data to 10 without the IC improves the description of the HERA charm production data at the expense of a worse fit to the inclusive DIS data and to the full data set. Including the nonperturbative charm contribution of the BHPS, SEA, or mixed type does not improve the fit to the EMC F 2c (x, Q), in contrast to the findings in [39]. It might be argued that a larger set of parametrization forms for the IC needs to be explored, as in the NNPDF method, to see if a better fit to the EMC F 2c (x, Q) could be reached. In the absence of control of experimental and (N)NLO theoretical systematic effects in the EMC F 2c data set, such an exercise again appears to be excessive. Indeed, when using a purely perturbative charm only, the NNPDF3.1 study [40] obtains a considerably worse χ 2 n /N pts = 4.8 for the EMC F 2c data set than our results quoted in Table I. After including a flexible "fitted charm" parametrization they arrive at a much better agreement with the EMC data sample, with χ 2 n /N pts = 0.93 and x c+c = 0.34±0.16% at Q c = m pole c = 1.51 GeV at 68% C.L. Their χ 2 /N pts values in Table 4.3 of [40] are somewhat better for the inclusive HERAI+II data set (1.16) and somewhat worse for the HERA charm SIDIS data set (1.42), compared to our 1.25 and 1.22 in Table I.
Some of these disparities are explained by non-identical PDF parametrization forms (positive-definite BHPS/SEA models in the case of CT14 IC, vs. the neural networks of NNPDF3.0), the general-mass schemes, and the choices of the mass parameters: Q c = m pole c =1.3, 1.275, and 1.51 GeV in the CT14, NNPDF3.0, and NNPDF3.1 studies, respectively. The preferred x c+c (Q) at Q = 1.51 GeV are smaller in the NNPDF3.1 framework than for CT14 IC in part because the evolved perturbative charm PDF is absent at this Q in NNPDF3.1. The S-ACOT-χ scheme that we use is at present the only ACOT scheme in which the massive coefficient functions are fully available to NNLO, or O(α 2 s ) [60]. NNPDF3.1 used a different mass scheme [37,39] and set to zero some O(α 2 s )/NNLO massive terms that are not available in that scheme [40]. We thus expect some differences between the schemes.
Another difference arises from the definitions of uncertainties. The current paper quotes 90% probability intervals obtained by scanning ∆χ 2 with respect to x , as explained in Sec. IV B. The NNPDF works quote their errors as symmetric standard deviations obtained from averaging over many replica fits, each of which is not a perfect fit and may deviate from the central fit by hundreds of units of χ 2 [113].
As an illustration, Fig. 14 compares the probability intervals on the momentum fractions from the CT14/CT14HERA2 and NNPDF3.0/NNPDF3.1 NNLO analyses. The left frame shows the CT14/CT14HERA2 90% probability intervals for x IC at Q = 1.3 GeV. The right frame shows the CT14/CT14HERA2 intervals for x c+c (Q) at Q = 1.51 GeV and superimposes the 68% C.L. uncertainties on the fitted charm (FC) copied from the NNPDF3.0 and 3.1 publications. Apart from the constant horizontal shift due to the Q c choice, without the EMC data, the CT14 and NNPDF probability intervals for x c+c are reasonably compatible, minding their non-equivalent definitions. [The upward shift in x c+c (Q) by ≈ 0.5% due to the choice of Q c , an auxiliary scale in a general-mass scheme, is of little physical significance, it is canceled up to O(α 3 s ) in the complete DIS cross section because of the compensating shift in ACOT subtraction terms.] Inclusion of precise LHC data sets helped to reduce the uncertainty in NNPDF3.1. The symmetric definition of the NNPDF3.1 errors allows a negative value of uncertain interpretation for x c+c at 68% C.L. if the EMC data are not included. A very small uncertainty on x c+c quoted by the NNPDF3.1+EMC fit is accompanied by the reduction in the global χ 2 by less than 13 units for 4300 data points when the EMC data are added into the fit, cf. Table 4.3 in Ref. [40]. Needless to say, the impact of the new experiments and assumptions on the uncertainty of x c+c warrants a further investigation.

CROSS SECTIONS AT THE LHC RUN II
Next, we will analyze the impact of the fitted/intrinsic charm (or the "IC", for short) on key observables at the LHC, assuming that the fitted charm does not strongly depend on the hard process at NNLO. [We argued in Section II that this assumption is not self-evident. We will nevertheless make it to investigate sensitivity of the LHC predictions.] Figure 15 illustrates dependence of the total cross sections for inclusive production of electroweak bosons W ± , Z 0 , and H (via gluon-gluon fusion) on the IC model and charm quark mass at the LHC √ s = 13 TeV. To provide a visual measure of the CT14NNLO uncertainty, each figure shows an error ellipse corresponding to CT14 NNLO at the 90% C.L. The W and Z inclusive cross sections (multiplied by branching ratios for the decay into one charged lepton flavor), are calculated by using the Vrap v0.9 program [16,17] at NNLO in QCD, with the renormalization and factorization (µ R and µ F ) scales set equal to the invariant mass of the vector boson. The Higgs boson cross sections via gluon-gluon fusion are calculated at NNLO in QCD by using the iHixs v1.3 program [114], in the heavy-quark effective theory (HQET) with finite top quark mass correction, and with the QCD scales set equal to the invariant mass of the Higgs boson. The first row of Fig. 15 shows predictions for W ± , Z 0 , and H 0 production cross sections in the five BHPS and SEA fits for m pole c = 1.3 GeV. Predictions for different values of the IC momentum fraction 0% < x IC < 3% and charm-quark mass 1.1 < m pole c < 1.5 GeV, obtained with the initial scale Q 0 = 1 GeV, are illustrated in the second and third rows of Fig. 15. The varied x IC values are indicated by the point color for each m pole c value. The central value predictions for the BHPS and SEA models are all within the CT14 NNLO uncertainties, with BHPS very close to the CT14 nominal fit. The impact of IC on these key LHC observables is mild. For BHPS, increasing x IC generally increases, and then reduces the W ± , Z 0 cross sections, and increases the Higgs cross sections. For SEA, increasing x IC reduces all cross sections.
The intrinsic charm may partially offset the variations in the electroweak cross sections due to the pole charm mass. As m pole c is increased from 1.1 to 1.5 GeV, the light-quark PDFs in CT14/CT14 HERA2 are mildly increased at x > 10 −3 and Q ∼ M Z , while the gluon is reduced at x > 0.1. As mentioned before, m pole c ≈ 1.5 GeV results in a worse fit to the CT14HERA2 data set, cf. the upper Fig. 6. For the LHC W/Z cross sections, increasing m pole c to 1.5 GeV results in two competing trends. On the one hand, 1.5 GeV leads to a somewhat better description of the total W and Z cross sections in Fig. 15, even though the changes are well within the CT14 uncertainty. This increase reflects larger u and d (anti)quark PDFs for m pole c = 1.5 GeV. On the other hand, the LHC data on high-p T Z-boson production [117][118][119] show contradictory preferences for the m c and x IC , depending on the collider energy [7 or 8 TeV] and the format of the data [absolute or normalized cross sections]. Our conclusion at the moment is that the LHC inclusive W and Z production cross sections may provide helpful correlated constraints on m c and x IC in the future. We may also consider more exclusive scattering processes [120][121][122][123][124][125][126][127][128] to look for evidence of the IC in the LHC environment. at NNLO and experimental points [115,116] are also shown.

VI. Z + CHARM-JET PRODUCTION IN PP COLLISIONS AT THE LHC
A suitable test scenario is given by the production of a Z boson in association with a charm jet, for which a CMS measurement at √ S = 8 TeV has been recently published in Ref. [129]. The corresponding calculation pp → γ/Z c is available at NLO in QCD, building on the important feature that the LO partonic process g + c → γ/Z + c (consisting ofŝ andt channel contributions) is directly sensitive to the initial-state charm distribution. Provided that the charm-quark transverse momentum is much larger than its mass, the NLO corrections to this process can be calculated working in the S-ACOT scheme [46,64,65]. Using this scheme enables one to neglect the charm mass throughout, while only making a small error of the order of 1/ ln M Z mc × m 2 c p 2 T [130]. The contributing subprocesses are given by gc → Zc (one-loop level production), q/g c → Zc q/g (light-flavour parton emission) (q = u, d, s) and gg → Zcc (charm pair production). 11 Another subprocess leading to charm-quark pairs in the final state is qq → Zcc. It is not regarded as a correction to gc → Zc, but it is an additional source of Zc events, and therefore taken into account at LO. There is one subtlety that concerns the Zcc final states. They are evaluated by retaining the charm-quark mass in order to regulate the gluon splitting singularity that would arise for massless collinear quarks. Taking all of these subprocesses, one then arrives at an NLO accurate description for the associated production of a Z boson and a single charm jet, as has been presented in Ref. [131] and implemented in the program MCFM [132]. To compare the impact of the different IC PDF fits, we use the MCFM calculation to generate the various Z+c-jet cross sections in the presence of intrinsic charm at NLO.
The main drawback of the fixed-order predictions is their limitation in describing effects that arise from multi-particle final states. One complication is due to the importance of jet production at higher orders, which enhances the size of the Z+charm-jet cross section especially for high-p T Z boson production. The inclusive cross section definition (Z+c-jet+X) employed by the CMS analysis makes it important to account for the contributions from more complex topologies like gq → Zccq or gluon-jet splitting to cc occurring at higher perturbative orders (i.e. in Monte Carlo physics language, later in the event evolution). The fixed-order approach will miss these multijet contributions, but we can invoke matrixelement plus parton-shower merging (MEPS) to study these effects. This can be particularly important if the final state is binned in a variable such as the Z boson transverse momentum, while a fixed (low) cut is placed on any jets in the event. We can also investigate at which point (in terms of the number of multileg MEs included), saturation (stabilization of the cross section) can be found.
Another complication stems from the fact that in an experimental environment, we are required to use a cross section definition, which is based on the detection of charm hadrons/objects in the event, i.e. charm tagging is involved one way or another to determine the inclusive Z+charm-jet rate. The theory-driven, parton-level definition employed in the fixed-order case cannot be applied here, as it ignores the evolution of the hard event to energy scales of the order of 1 GeV, where the measurement takes place. In this context and, especially, for the identification of specific particles/objects -as in our case, charm jets -aspects of multi-particle production (beyond hard jets) therefore need to be taken into account to arrive at a more realistic simulation. For our studies, we will rely on the parton shower to describe the fragmentation of the charm partons [133], and assuming factorization of the initial-state and final-state QCD radiations as a reasonable approximation. The cross section based on charm tagging will be affected by parton showering. Thus, we have to deal with contributions emerging from Z+non-c partonic processes because the g → cc splittings have the potential of turning a Z plus light-flavour jet into a Z plus c-jet contribution. This additional source of Z+charm events enhances the size of the measured cross section. However, this enhancement simply serves to dilute the impact of any intrinsic charm, since in most cases it emerges from initial states not involving a charm quark, i.e. the enhancement comes from final-state gluon splitting into a cc pair. The rate for this enhancement depends on both the charm-jet transverse momentum threshold and the number of jets in the final state.
For these reasons, one cannot ignore the multi-particle aspects when dealing with realistic scenarios. We therefore generate predictions using the LO matrix-element plus partonshower merging (MEPS@LO) approach [134], adding additional jets and subsequent parton showering, and requiring the presence of a charm jet in the final state. The MEPS@LO approach allows us to estimate the impact of the higher-order radiative corrections and charm tagging at the same time. Using the various IC models, we can examine (on a quantitative level) to what extent the multi-particle effects alter the outcome of the NLO calculations provided by MCFM. All MEPS@LO predictions presented here have been obtained from the Monte Carlo event generator Sherpa-2.2.1 [135]. To perform the charm tagging in the Sherpa simulations, we rely on the flavorful version of the anti-k t jet algorithm as implemented in FastJet [136]. We generate Z+jets samples in the five-flavour scheme (massless c and b quarks) involving tree-level matrix elements for Z+0, 1 parton up to those for Z+n ME partons where n ME denotes the maximum outgoing-parton multiplicity of these matrix elements. Three Sherpa samples are provided, namely for n ME = 1, 2, 3, using a merging cut of Q cut = 20 GeV. Each Sherpa Nj prediction is then drawn from the respective Z+jets sample with n ME = N.
The simplest observable to look at is the inclusive Z+charm-jet cross section. Hence, we start by presenting a summary of cross section predictions in Table II. In both types of calculations (fixed order and MEPS@LO), we employ kinematic requirements that are similar to those utilized by the CMS analysis [129]. 12 Most notably, we impose the following  kinematic requirements on the two leptons from the Z boson decay: p T,ℓ > 20 GeV, |η ℓ | < 2.1 and 71 GeV < m ℓℓ < 111 GeV. Jets are defined by using the anti-k t algorithm with a size parameter of 0.5 and threshold requirements reading p T,jet > 25 GeV as well as |η jet | < 2.5. Table II shows that the predictions from the two standard PDFs (CT14NNLO and NNPDF3.0) agree very well. All CT14 IC models lead to an increase of the Z+charm-jet cross section varying from about 2% for the specific choice of using BHPS3 to almost 20% for the SEA2 model. On the contrary, the fitted charm PDFs of the NNPDF group [39,112] lead to a small reduction of the total cross section, however by no more than 5%. The results from Table II also confirm the rise of cross sections owing to the inclusion of multijet contributions. This increase can grow as large as 10%. From a fixed-order point of view, the Sherpa 1j calculation is LO-like while the Sherpa 2j calculation is closest to the one provided by MCFM. The largest differences with respect to MCFM lie in Sherpa's neglect of virtual contributions that are non-Sudakov like and the usage of a dynamical plus local scale setting prescription. 13 The Sherpa 3j computation then goes beyond the MCFM (cf. Ref. [129]), we refrain from comparing our results directly to the experimental data for reasons such as unapplied/unknown hadronization corrections and neglecting certain ∆R constraints. 13 We note that the Sherpa 2j calculations can be made even more MCFM-like by relying on Sudakov reweighting but applying no parton showers at all. These modified Sherpa predictions show good agree- calculation, resulting in an additional but smaller increase with respect to the Sherpa 2j cross sections. In other words, we observe the expected saturation effect that stabilizes the Z+c-jet rate with increasing n ME . As in the fixed-order case, the CT14 IC models enhance the Sherpa cross sections by different amounts. For a specific model, the predicted gains are of similar size among the different Sherpa Nj calculations (as indicated by the numbers in parentheses in Table II), but turn out to be smaller when compared to the respective fixed-order result. The MEPS@LO predictions therefore show the expected dilution of the IC signals as previously described. Furthermore, we can take this as evidence for similar mitigation effects applying to experimental signatures for intrinsic charm. The total inclusive cross section as measured by CMS, σ(pp → Zc + X) × BR(Z → ℓ + ℓ − ) = 8.6 ± 0.5 (stat.) ± 0.7 (syst.) pb, comes with an overall relative uncertainty of 10%. ment with the cross sections predicted by MCFM though they are still larger by about 2%. This cross section is larger than any of the predictions in Table II. With this rather large value, we cannot yet draw any conclusion regarding the preference or exclusion of the various IC models. For example, if we assume that the baseline CT14 prediction describes the data, the SEA2 model, which predicts the largest relative cross section change among all IC models, would only occur at the upper edge of the allowed (2σ) range (neglecting the impact of PDF and theory uncertainties for a moment). However, the various intrinsic charm models affect the low and high x regions differently, making it worthwhile to investigate the effects on differential cross sections as well. As mentioned earlier, the transverse momentum distribution of the Z boson in association with a charm jet is a suitable candidate because larger x values predominantly affect the high p T region. Focusing on different p T regions may therefore increase our chances to distinguish certain IC models from each other. Figures 16 and 17 show MCFM predictions of the differential Z boson p T cross sections at the LHC, for energies of 8 TeV and 13 TeV, respectively. Apart from presenting the p Z T distributions themselves, we also depict the respective ratios taken with respect to the CT14NNLO result. We furthermore use the panels on the right in Figures 16 and 17 to present a similar set of plots obtained by using various PDFs from the NNPDF group, namely their current default version, NNPDF3.0, also serving as the reference curve in the lower part of the right panes, and their associated fitted charm PDFs with and without accounting for the EMC data. These NNPDF plots also contain the CT14 baseline predictions (including their PDF uncertainties) to allow for direct comparison between both PDF families. show the ratios between the Sherpa 3j prediction using CT14NNLO and those using the IC models.
The results of Figures 16 and 17 reveal the existence of sizable deviations between the predictions from the standard PDFs and the IC models (for both families). The BHPS intrinsic charm fits produce larger cross sections for high Z transverse momenta, while the PDFs using the SEA parametrization affect the cross sections fairly equally at all values of p Z T , and in a similar way at both 8 TeV and 13 TeV predictions. In particular, the SEA2 fit yields increases of the order of 20%. Regarding the BHPS models, the critical issue is the reach of the LHC data into regions of higher x (corresponding to large values of p Z T ) where the enhancement in the BHPS models becomes significant. At 8 TeV, the effects can be up to 100% higher than the baseline; they however occur in a region without data (for p Z T > 500 GeV). At 13 TeV, we deal with smaller x values on average and therefore observe smaller deviations (dropping by nearly a factor 2) for the corresponding BHPS predictions. We also note that the relative changes predicted by the fitted charm PDFs of the NNPDF group resemble those of the BHPS fits for the CT family. This resemblance is found at both collider energies, for which we also observe good agreement between the central predictions of NNPDF3.0 and CT14NNLO.
As discussed previously, we expect the sensitivity to the intrinsic charm component to de- crease in a realistic multijet environment. The p Z T distributions provided by the MEPS@LO method for the various PDFs are presented in Figures 18 and 19, to be compared with Figures 16 and 17 depicting the corresponding MCFM results. To support a direct comparison, the main panels of Figure 18 also contain the MCFM prediction for CT14NNLO. While there are no large deviations between the Sherpa and MCFM predictions for lower p Z T values, the Sherpa predictions show the expected hardening in the tail of the p Z T distributions. In the MEPS@LO simulations, the IC increases the cross sections in the same way as in the fixed-order case, although by a smaller factor (roughly half as much), which is most prominently visible in the associated ratio plots.
Apart from reconfirming the dilution effect, Figure 19 provides us with additional information. First, the Sudakov (low p Z T ) region is described in a more sophisticated and therefore robust way (as a result of the inclusion of resummation effects). Second, regardless of whether CT14NNLO (Figure 19-left) or BHPS2 is used as reference (Figure 19-right), the inclusion of additional layers of multileg matrix elements leads to relative enhancement and  saturation effects of similar size at larger p Z T values. This is an expression of the fact that although the intrinsic charm models investigated here do change the initial conditions of the charm content in the proton, they do not alter the nominal QCD evolution. The parton shower evolves in the same way as encoded by the DGLAP theory in the absence of any intrinsic charm.
Similarly to the case of the Z+c-jet cross section, the CMS data for the p Z T distribution [129] can be used to estimate the current potential for discriminating possible intrinsic charm models. The CMS measurement provides cross sections for three different p Z T /GeV bins, which are shown in the upper part of Table III, together with their associated relative uncertainties. These uncertainties are to be compared with the size of the deviations induced by the intrinsic charm fits with respect to the CT14 baseline. According to Figures 16 and  18, we can focus on the BHPS2 and SEA2 predictions, as only those feature differential rates significantly exceeding the uncertainty range of the CT14 prediction. However, as shown in the lower part of Table III, the deviations generated by both the BHPS2 as well as the SEA2 model do not exceed the 1σ variation of the data, in particular if the dilution effect is taken into account as simulated by MEPS@LO. Thus, none of these changes reach a magnitude that is distinguishable from the experimental and theoretical systematic errors at 8 TeV. The discriminating power of the current CMS data is simply not sufficient to test the IC models, either in terms of the differential p Z T cross section or in terms of the total Z+charm-jet cross section. 14 Owing to the rather low impact of current LHC data, it is important to better understand the prospects for new measurements of detecting or excluding a high-x IC component. To this end we extrapolate what we have learned at 8 TeV to the case of the 13 TeV LHC. The CMS result for 19.7 fb −1 of data at 8 TeV extends up to a Z boson transverse momentum range of 200 GeV. The last bin is fairly wide, from 60 GeV to 200 GeV, and its associated differential cross section has been measured as ∆σ/∆p Z T = (0.017 ± 0.003) pb/GeV, i.e. it is reasonable to assume that cross sections as low as 0.01 -0.02 pb/GeV can be measured with ∼ 20 fb −1 of integrated luminosity. In Figure 16, a cross section of this size corresponds to a p Z T value of about 120 GeV, which translates into x ∼ 0.03 on average. Thus, current measurements probe relatively low values of x compared to the range (x ≥ 0.1 -0.2) where the BHPS models start to have a significant impact (as shown in Figure 8). The cross section for Z+charm production of course is larger at 13 TeV, but for the same low cross section target of 0.01 pb/GeV, the accessible p Z T range would only be extended by 30 GeV (according to Figure 17) pulling the mean x towards 0.02, which means we would not even achieve the same sensitivity as for the 8 TeV case. To reach a similar x range would require a Z transverse momentum of the order of 200 GeV corresponding to a cross section of about 0.002 pb/GeV. One therefore needs an integrated luminosity of about 100 fb −1 , in order to determine this cross section with an accuracy comparable to the 8 TeV case. In other words, it will take the full Run 2 cycle to barely get a first 2σ sign of deviations at p Z T ∼ 200 GeV or probe transverse momenta of the order of 300 GeV. Needless to say that definitive confirmation/exclusion will require us to go considerably beyond the Run 2 luminosity budget.
The challenging environment for the Z + c analysis forces us to search for ways to increase the impact of an intrinsic charm component on the Z+c-jet cross section. As this cross section is diluted by the presence of the radiative corrections, for example, limiting the number of jets in the event could reduce this dilution. The Z+c-jet rate could also be measured as a function of the leading (charm) jet transverse momentum, which in fact has been carried out by CMS in the same publication. Our studies suggest that this differential cross section is somewhat more sensitive to the intrinsic charm modeling investigated here, but its sensitivity must be weighed against the size of the relative uncertainties on the measurement of the charm jet p T , in a similar fashion as shown for p Z T in Table III. In addition, deviations are also found for Z boson rapidities outside the central phase-space region, such as might be measured at LHCb [128].

VII. SUMMARY AND CONCLUSIONS
We have explored the possibility of having a sizable nonperturbative contribution to charm parton distribution function (PDF), i.e., the intrinsic charm (IC) quark component, in the proton, using the CTEQ-TEA (CT) global analysis. In Sec. II, we reviewed the theoretical framework used in the CT global analysis, and discussed the conditions under which our formalism, Eq. (1), can better approximate the QCD factorization theorem, Eq. (2).
The notion of "intrinsic charm" refers to contributions to charm quark production and scattering that arise besides twist-2 "perturbative" contributions. In DIS, the twist-4 cross sections for charm production may numerically compete with "perturbative" twist-2 cross sections at a high enough order in α s . For example, in Fig. 3, we show the relevant squared amplitudes for a DIS structure function F (x, Q) from the γ * + gg → c +c process. The flavor-creation diagrams F (n) h,gg render most of the twist-4 charm production rate in the HERA kinematical region (Q m c ). But, at very high photon virtualities, Q 2 ≫ m 2 c , their dominant part is approximated in a variable-flavor number scheme by a twist-2 coefficient function c (k) h,h convoluted with a universal charm PDF c(x, Q). A non-zero boundary condition for c(x, Q) at Q = Q c ∼ m c is obtained by perturbative matching from light-parton nonperturbative twist-two and twist-four functions, such as f g/p (x, Q c ) and f gg/p (x 1 , x 2 , Q c ).
In the context of the phenomenological PDF analyses, on the other hand, the "intrinsic charm" PDF is often conflated with a "fitted charm" PDF parametrization that plays a dual role of the approximant for the above power-suppressed contribution to charm scattering and of a parametric surrogate for unrelated radiative contributions that were not explicitly included. At the moment the fitted PDF is determined solely using the fixed-order convolutions with the twist-2 coefficient functions, without including explicit twist-4 terms. While the "fitted charm" PDF provides a good description of the cumulative QCD data in the CT fit, care is necessary when making predictions for new processes based on its parametrization, as it may absorb a host of process-dependent corrections, notably the contribution of DIS-specific twist-4 coefficient functions like in Fig. 3. We, as well as the other global analysis groups, treat the "fitted charm PDF" obtained this way as though it is mostly process-independent, until it is demonstrated otherwise.
For example, in neutral-current DIS charm production the twist-4 charm cross section is of the same order in the QCD coupling strength as the NNLO twist-2 one. To estimate the magnitude of the twist-4 IC cross section from the DIS data, using its model given by the fitted charm, the twist-2 DIS contributions in the fit must be evaluated at least to NNLO. Furthermore, it is necessary to study the contributions from the strange (and bottom) PDF, dependence on the charm quark mass (m c ), and to accurately implement suppression of charm production at the mass threshold. In the case when low-Q fixed-target data are included, the IC component must be further discriminated from the 1/Q 2 and nuclear-target effects.
Hence, in this study, we have used both the CT14 NNLO and CT14HERA2 NNLO analyses, differing mainly in their strange PDFs. CT14HERA2 has a softer strange quark component than CT14 at most x values. We have carried out a series of fits with a varied charm quark pole mass m c between 1.1 and 1.5 GeV, within the preferred m pole c range of our global fits, see Fig. 6.
The NNLO heavy-quark mass effects are implemented in our calculation using the S-ACOT-χ factorization scheme. 15 In Sec. II, we have given detailed arguments showing that it is a self-consistent and sufficient scheme for predicting massive-quark DIS contributions both in the twist-2 and twist-4 channels.
The charm content in a hadronic bound state, quantified by an operator matrix element identified with the charm PDF, can in principle be predicted by QCD. We examine which "intrinsic charm" models predict the fitted charm PDF compatible with the global QCD data. Two generic types of the charm models introduced in Sec. III, a valence-like BHPS model and a sea-like SEA model, predict a non-zero c(x, Q 0 ) at large x and across all x, respectively. The BHPS model is solved either approximately in the BHPS1 and BHPS2 PDF sets, or exactly in the BHPS3 set. To better predict the PDF ratios of charm to up and down PDFs, in the BHPS3 model we also allowed for small intrinsic contributions to theū andd (anti-)quarks generated from the |uuduū and |uuddd Fock states, included together with the charm intrinsic contribution. Though we did not present its details, we have also studied a mixed model of BHPS and SEA and arrived at similar conclusions. Figure 5 shows that, at Q 0 = 1.3 GeV, the charm quark momentum fraction x IC , as defined in Eq. 13, is found to be less than about 2% and 1.6%, for the BHPS IC and SEA IC models, respectively, in the CT14NNLO analysis, at the 90% C.L. We note that by its definition, x IC is evaluated at the initial scale Q 0 . It is to be distinguished from the full charm momentum fraction x c+c at Q > Q c , which rapidly increases with Q because c(x, Q)+c(x, Q) also includes the perturbative contribution. The dependence of the outcomes on m pole c was reviewed in Sec. IV C, and the resulting BHPS and SEA PDFs and parton luminosities, as well as Q dependence of x c+c , were explored in Section IV D.
A significant IC component in the proton wave function could influence observables measured at the LHC, either directly through enhanced cross sections via the charm PDF, or indirectly via the momentum sum rule leading to a change in the momentum fraction carried by the gluons. Modifications in the light-flavor PDFs are generally mild in the considered BHPS/SEA models, although the gluon-gluon luminosities can be suppressed at the highest final-state invariant masses M X , as observed in Fig. 10. The allowed momentum fraction x IC is correlated with the charm pole mass m pole c , especially in the SEA model. When the charm PDF is purely perturbative, the inclusive Z cross section increases as m pole c increases, due to the largerū andd PDFs that compensate for the smaller perturbative charm PDF 15 The massive NC DIS perturbative coefficients are known in their entirety to NNLO in the S-ACOT-χ [60], TR' [48], and FONLL-C [49] schemes. In contrast, some of these NNLO coefficients are still unknown in the "fully massive" ACOT scheme [46] and its FONLL equivalent [37,38] adopted in NNPDF3. 1. contribution. We also observe reduction in g(x, Q) at large x, and consequently some reduction in cross sections sensitive to large-x gluon scattering. For example, increasing m pole c from the nominal 1.3 to 1.5 GeV increases the W/Z inclusive total cross sections at 13 TeV, reduces the normalized high-p T Z production cross section at the LHC 7 TeV, and has vanishing effect on the gg → H 0 cross sections, see Sec. V. These changes can be partly offset by introducing the IC, possibly at the expense of some tension with the non-LHC fitted experiments, and generally within the regular CT14 PDF uncertainty.
There is much discussion in the literature about the impact of the EMC measurement [41] of semi-inclusive DIS charm production on the intrinsic charm PDF. Although our standard analysis does not include the EMC data, we have examined their impact in several IC models. Section IV F argues that fitting the EMC data is not expedient, their persistent tension with the other fitted data sets may reflect the systematic errors that were not documented in the EMC publication. The level of (dis)agreement with the purely perturbative charm and the exclusion limits on the intrinsic charm depend on the assumed magnitude of systematic effects in the EMC measurement. As shown in Table I, even without the IC contribution, the χ 2 /N pts of the EMC data varies from about 3.5 to 2.3 when it is excluded or included with a large statistical weight in the CT14 fits. Including the intrinsic charm component does not significantly change χ 2 /N pts for the EMC. For the BHPS models, including the EMC data with the nominal errors reduces the tolerated range of x IC by about a factor of two. The impact of EMC data is small within the SEA model.
Besides the LHC electroweak boson production cross sections, we examined the implications of the IC for associate production of Z boson and charm-jet at the LHC, and summarized our findings in Table II and Figs. [16][17][18][19]. A fixed-order calculation for Z + c production, MCFM at NLO, was compared to a merged parton showering calculation in Sherpa, which also generates charm jets in the final state via gluon splittings. In general, in a fixed-order calculation for Z + c, the various IC models predict enhanced rate in the transverse momentum distribution of a Z boson (p Z T ) [128]. The SEA models tend to predict a higher differential cross section across all p Z T , while the BHPS models suggest the increased rate only at the highest p Z T . The predictions based on the NNPDF3IC and NNPDF3IC (no EMC) PDFs are close to our BHPS3 and BHPS2 predictions, respectively, they predict a larger rate in the high p Z T region. Inclusion of the final-state parton showering typically dampens the fixed-order enhancement induced by the IC contribution, as can be observed from the comparison of Sherpa to MCFM predictions. The dampening is mainly attributed to the gluon-splitting contributions in the final state which reduce the relative impact of the IC contribution in the hard p Z T tail, especially for the predictions from the BHPS models.
The analysis of QCD factorization indicates that the power-suppressed "intrinsic" component in semi-inclusive DIS charm production may be comparable in magnitude to some NNLO and N3LO leading-power contributions. Hence, a serious study needs to be carried out at least at the NNLO, such as in this work. (It is not possible to draw a definite conclusion from an NLO analysis.) As of today, the experimental confirmation of the IC component in the proton is still missing, and data from far more sensitive measurements are required. An analysis of very low-Q fixed-target data, such as the one presented at NLO in Refs. [33,34], must meet the challenge of the reliable separation of the IC from the other relevant factors, including higher-order twist-2 contributions, the 1/Q 2 terms, m c dependence, and nuclear effects. The constraints on the IC from the higher-energy data are largely compatible between the CT14 IC and NNPDF3.x analyses [39,40]. Our limits on x c+c are moderately more conservative than those of NNPDF3.1, as we do not include the EMC F 2c data and acknowledge competing preferences for m pole c and x c+c among the various non-LHC and LHC experiments, as outlined in Secs. IV C, IV E, and V. Ultimately, a combination of high-luminosity measurements at the Large Hadron Collider, such as Z + c production, and charm SIDIS at the Electron-Ion Collider [137] will be desirable to test intrinsic charm scattering contributions at NNLO and beyond.