#
CT14 Intrinsic Charm Parton Distribution Functions

from CTEQ-TEA Global Analysis

###### Abstract

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.

###### pacs:

12.15.Ji, 12.38 Cy, 13.85.Qk## 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 Dulat et al. (2016) 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 , for the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) evolution equations Gribov and Lipatov (1972a, b); Lipatov (1975); Dokshitzer (1977); Altarelli and Parisi (1977), 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 Moch et al. (2004); Vogt et al. (2004); Sanchez Guillen et al. (1991); van Neerven and Zijlstra (1991); Zijlstra and van Neerven (1991, 1992); Laenen et al. (1993); Riemersma et al. (1995); Buza et al. (1996); Anastasiou et al. (2003, 2004). Next-to-leading order (NLO) is used only for inclusive jet data from the Tevatron and the LHC and for deep-inelastic scattering (DIS) charged-current data from HERA and fixed-target experiments. The NNLO predictions for these processes Berger et al. (2016); Currie et al. (2017a, b) were not available or incomplete at the time of the CT14 study, and we have argued Gao et al. (2014); Dulat et al. (2016) 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 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 of order of the charm mass. In the CT14 analysis, the charm quark and antiquark PDFs were turned on at the scale GeV, with an initial distribution consistent with NNLO matching Buza et al. (1996, 1998) to the three-flavor result. At higher , 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 and perturbatively evolved to the experimental data scale .

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 Fock state of the proton wave function Brodsky et al. (1980, 1981); Pumplin (2006); Chang and Peng (2011); Blümlein (2016); Brodsky et al. (2015) and predicted by meson-baryon models Navarra et al. (1996); Paiva et al. (1998); Steffens et al. (1999); Hobbs et al. (2014). On the other hand, Refs. Jimenez-Delgado et al. (2015, 2016) 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. Brodsky and Gardner (2016). This is to be compared to the earlier CT10IC study Dulat et al. (2014), 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 Ball et al. (2016a, 2015a, b, 2017) established a smaller fitted charm momentum fraction. NNPDF determined a fitted charm momentum fraction equal to )% at 68% C.L. just above the charm mass threshold, with the charm quark pole mass taken to be 1.51 GeV Ball et al. (2017), and equal to )% when the EMC data Aubert et al. (1983) on SIDIS charm production were included.

The current paper revisits the issue in the context of the CT14 analysis Dulat et al. (2016), also including more recent advances that were made in the follow-up CT14HERA2 study Hou et al. (2017a). It updates the previous work Dulat et al. (2014) on fitting the charm PDFs based on the CT10 NNLO framework Gao et al. (2014), as well as the CTEQ6.6 IC study Pumplin et al. (2007) 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 and 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 Jimenez-Delgado et al. (2015, 2016); Lyonnet et al. (2015); Ball et al. (2016b, 2017). 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 , 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 Aivazis et al. (1994a, b); Buza et al. (1998); Thorne and Roberts (1998); Martin et al. (2010); Forte et al. (2010), the order of the calculation, the assumed charm mass , 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 in the absence of the nonperturbative charm has been addressed at NNLO in the CT10 NNLO framework Gao et al. (2013) and in other references Alekhin et al. (2013); Harland-Lang et al. (2016); Ball et al. (2016b); Bertone et al. (2016); Alekhin et al. (2017); Gizhko et al. (2017). In the context of the CT10 analysis Gao et al. (2013), the general dependence on the charm quark mass was studied, and a preferred value of GeV was obtained at 68% C.L., where the error is a sum in quadrature of PDF and theoretical uncertainties. Here, denotes the running mass of the charm quark, defined in the modified minimal-subtraction () scheme and evaluated at the scale of . This value, constrained primarily by a combination of inclusive and charm production measurements in HERA deep-inelastic scattering, translates into the pole mass GeV and GeV when using the conversion formula in Eq. (17) of Ref. Chetyrkin et al. (2000) 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 Bigi et al. (1994); Beneke and Braun (1994); Beneke (1999). In this light both converted values are compatible with the value of 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 Dulat et al. (2016) and CT14HERA2 sets Hou et al. (2017a) 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 and . Hence, accurate predictions of the and parton distributions will be relevant to various important LHC measurements, such as production of , , and Higgs boson, or associated production of a charm jet and a .

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 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.3. In Sec. IV.6, 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 , and Higgs bosons at the LHC in Sec. V, and on a charm jet production associated with a boson at the LHC in Sec. VI. Finally, our conclusions are presented in Sec. VII.

## Ii QCD factorization with power-suppressed charm contributions

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
short-distance 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
( ), heavy quarks ( and ), and gluons
.^{1}^{1}1Without loss of generality, we focus on a situation when
neither top quarks nor photons are classified as nucleon’s partonic
constituents. A parton 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
, or for short, depends on the parton’s flavor
and, ultimately, the parton’s mass.

A charm quark with mass 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 at an initial energy scale of order 1 GeV and evolve the PDFs to higher energies using the DGLAP equations Gribov and Lipatov (1972a, b); Lipatov (1975); Dokshitzer (1977); Altarelli and Parisi (1977). For the charm and anticharm contributions, on the other hand, the respective PDFs at such low are not mandatory. Only some QCD factorization schemes introduce them, with the goal to improve perturbative convergence at scales much larger than . 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 with 3 and 4 active flavors, at a relatively low momentum transfer that is comparable to the mass of the charm quark. Our considerations can be extended readily to situations with more than four active flavors, and to higher values. Moreover, among the experimental processes included in the global QCD analysis, the neutral-current DIS is the most sensitive to charm scattering dynamics Gao et al. (2013); Harland-Lang et al. (2016); Ball et al. (2016b); Bertone et al. (2016); Alekhin et al. (2017); Gizhko et al. (2017) with the rest of the processes providing weaker constraints. Therefore, it is natural to focus on DIS as the starting point.

### ii.1 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:

(1) | |||||

This is a standard convolution formula, consisting of the coefficient
function
and the PDFs dependent on the light-cone
partonic momentum fraction and factorization scale of
order (set to coincide with the renormalization scale to simplify
the notation). The index denotes the initial-state parton’s flavor, running
from , corresponding to the gluon, to the number of
active quark flavors assumed in the QCD coupling
strength and the PDFs Implicitly,
summation over quarks and antiquarks is assumed. We reserve the index
“ for a heavy-quark flavor, in DIS charm production.^{2}^{2}2
Beyond the NNLO accuracy considered in this paper, DIS includes
contributions with both and quarks. Treatment of such
contributions in the ACOT formalism is explained in Refs. Guzzi et al. (2012); Wang (2015). The superscripts in both and
emphasize that their perturbative coefficients are computed up to a fixed order
of .

Let us highlight several aspects of this formula. First, , the number of active flavors, is not measurable, it is a theoretical parameter of the renormalization and factorization schemes chosen for the perturbative calculation. should be distinguished from Tung et al. (2007); Guzzi et al. (2012), the number of (anti-)quark species that can be physically produced in the final state in DIS at given collision energy. The optimal value of is chosen as a part of the QCD factorization scheme to optimize perturbative convergence. 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 in the S-ACOT- scheme Aivazis et al. (1994b); Collins (1998); Kramer et al. (2000); Tung et al. (2002), a general-purpose factorization scheme for lepton-hadron and hadron-hadron scattering processes. For neutral-current DIS, were derived in this scheme up to , or NNLO Guzzi et al. (2012). Figure 1 is reproduced here from Ref. Guzzi et al. (2012) and shows the Feynman diagrams and notations for the perturbative coefficients of the “charm production” structure function 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 and not as important numerically, can also be found in Ref. Guzzi et al. (2012).

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,
The input charm mass can be either the
mass , or the pole mass
: the two are related by NNLO perturbative
relations Chetyrkin et al. (1998, 2000), both are
implemented in CT14 PDFs.^{3}^{3}3
The past CTEQ-TEA analyses traditionally used as an input, but may be preferable in future precise calculations.
The pole mass cannot be used to arbitrarily high accuracy because
of nonperturbative infrared effects in QCD, related to
the fact that the full quark propagator has no pole
because of the quark confinement Patrignani et al. (2016).
As
are evolved upward from the initial scale , they
are converted from to , and from to at the
corresponding switching points . The perturbative
coefficients of are converted
concurrently to preserve the factorization
scheme invariance at each order of . The CT14 analysis
switches from to exactly at the heavy quark mass;
so for the charm quark the switching takes place at the energy scale .

In this conventional setup, we assume a zero charm PDF, , for at the initial scale slightly below , and obtain a small non-zero for at scale via perturbative matching. Of course, is arbitrary, we could equally choose a value below and then expect a non-zero charm PDF also at . This alternative suggests the possibility of including a non-zero initial charm PDF parametrization, or the “fitted charm” parametrization, at the initial scale that would now correspond to . However, if the charm quarks are produced exclusively from perturbative gluon splittings, the dependence on the fitted cancels up to the higher 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 with massive quarks that follows from the QCD factorization theorem Collins (1998, 2013):

(2) |

Eq. (2) underlies all modern computations for the inclusive DIS observables, in the factorization schemes with fixed or varied values. The convolution of with in Eq. (2) includes all “leading-power” radiative contributions that do not vanish when the physical scales , 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 is approximated just up to order

This means that, in the all-order factorization theorem (2), , 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 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 leading-power perturbative contributions beyond the order . The , or NNLO, radiative contribution to neutral-current DIS heavy-quark production is large numerically. If a global fit is done at NLO, as in Refs. Pumplin et al. (2007); Jimenez-Delgado et al. (2015); Ball et al. (2016b), it prefers an augmented fitted charm 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 right-hand 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 with integer (“twist-4”, see, e.g., Jaffe (1983); Jaffe and Soldate (1982)). In the DIS scattering of charm quarks, the lowest power-suppressed contribution also includes terms of order Brodsky et al. (1980, 1981); Collins (1998). The latter term clearly does not vanish with increasing and, furthermore, at very high it is enhanced logarithmically and behaves as with due to contributions from collinear scattering. The power-suppressed charm contribution, once introduced at low scale , will survive to the much higher scales relevant to the LHC.

### ii.2 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 leading-power, or twist-2, terms Guzzi et al. (2012). It is useful to compare the relevant Feynman graphs in the factorization scheme, the most appropriate scheme to use in the threshold kinematical region, where is comparable to , and in the scheme, which is most appropriate at , where the charm density has the most physical interpretation.

First, recall that in the scheme all subgraphs containing heavy-quark propagators are assigned to the Wilson coefficients and not to the PDFs . 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 . The explanation for this is that the 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 and . The nonperturbative subgraphs with virtualities much less than contain only light-parton propagators, as those are renormalized in the scheme.

A twist-4 hard-scattering matrix element for can be thought of as a twist-2 hard-scattering matrix element connected to the parent hadron by an additional light-parton 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 . Insertion of two QCD vertices suppresses the twist-4 cross section by a power of 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 with of order .

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 Braun et al. (2010); Ji and Belitsky (2015). The complete analysis of twist-4 is lengthy – we will refer to the vast literature on the subject, including Refs. Gross and Treiman (1971); Anikin and Zavyalov (1978); Jaffe (1983); Jaffe and Soldate (1982); Ellis et al. (1983, 1982); Balitsky et al. (1990); Qiu and Sterman (1991a, b); Jaffe and Ji (1992); Jaffe (1996); Müller et al. (1994); Blümlein et al. (1999); Geyer et al. (1999); Geyer and Lazar (2001).

We further note that, in the limit 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 containing the ladder contributions. denotes a two-particle irreducible (2PI) part (in the vertical channel) of the structure function . 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 Curci et al. (1980). It is composed of 2PI subgraphs and (without an upper index), where and are coupled to the virtual photon and target hadron, respectively. The decomposition in terms of and for twist-2 also appears in the Collins’ proof of QCD factorization for DIS with massive quarks Collins (1998).

The ladder graphs are different in the and schemes. The scheme introduces additional terms with heavy quarks that approximate the leading contribution in the limit. In Fig. 1, these ladders correspond to the contributions proportional to the “flavor-excitation” Wilson coefficient functions . Such terms are absent in the scheme, and their purpose is to resum collinear logs from higher orders with the help of DGLAP equations. In this case both the light- and heavy-parton subgraphs are renormalized in the scheme. Importantly, apart from a finite renormalization of , the perturbative expansions of the structure functions in the and schemes are equal up to the first unknown order in – 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 [with an upper index “(4)”], connected to the twist-2 kernels in the upper part via a “mixed-twist” kernel containing a real heavy-quark emission. As is connected to by four propagators, at it scales as . Since it includes loop integrals with massive quark propagators , the momentum scale can be either or ; but the term is less suppressed than . [It is crucial that two large QCD scales, and , are present, in contrast to the massless-quark case.] On the other hand, apart from the replacement of by , the second ladder has the structure of the first one.

#### ii.2.1 Factorization for twist-2 contributions

We assume that the Feynman diagrams in Fig. 2 are unrenormalized and indicate this by a subscript “0”. Ref. Collins (1998) shows how to recast the full sum of twist-2 diagrams into a factorized convolution

(3) |

by recursively applying a factorization operator and renormalizing the UV singularities. 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 operator is to replace the exact ladder graph by a simpler, factorized expression which provides a good approximation to the full graph in the limit, and which is valid up to a power-suppressed remainder . In particular, 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 operators to all orders, one demonstrates factorization for in either the scheme or the scheme of the Aivasis-Collins-Olness-Tung (ACOT) class Aivazis et al. (1994b). By its construction, the remainder is of order

(4) |

with or .

While the operation in the 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 , may include additional powers of that vary among the conventions Collins (1998); Guzzi et al. (2012). The ambiguity in gives rise to several versions of the ACOT-like schemes, all equivalent up to a higher order in . The form of may be even made dependent on the type and order of the scattering contribution: some choices for , such as the one made in the SACOT- scheme Kramer et al. (2000); Tung et al. (2002); Guzzi et al. (2012), simplify perturbative coefficients and enable fast perturbative convergence.

In a practical calculation of a twist-2 cross section illustrated
by Fig. 1, the operation defines the prescription
for constructing the perturbative coefficients
of Wilson coefficient functions from the structure functions
computed in DIS on a partonic target
. Here denotes an (anti)quark struck by the virtual photon.^{4}^{4}4Up to NNLO, we use a simplified decomposition of the neutral-current
DIS structure function over the quark flavors probed by the virtual
photons: ,
where is the (anti)quark’s electric charge Guzzi et al. (2012). The decomposition of the ACOT structure functions for higher orders was derived in Ref. Wang (2015).
The parton-scattering structure functions, coefficient functions,
and PDFs are expanded as a series in :

(5) |

where are perturbative coefficients Buza et al. (1996) of operator matrix elements for finding a parton in a parton A perturbative coefficient of the Wilson coefficient function at order can be found by comparing the perturbative coefficients on the left and right sides of

(6) |

The comparison does not specify the form of the perturbative coefficients with an initial-state heavy quark; those are specified by at each order and re-used in exactly the same form in all occurrences of in the contributions of orders and higher. The freedom in selecting affects and not the partonic PDF coefficients that remain defined in the scheme. With such self-consistent definition, the dependence on cancels up to the first unknown order in , as it was verified numerically up to in Ref. Guzzi et al. (2012).

#### ii.2.2 Factorization of twist-4 contributions: a sketch

Going back to Fig. 2, we recall that, while in the twist-2 factorization formula (3) the subgraph is counted as a part of the remainder , diagrammatically, it is attached to the upper ladder subgraphs in exactly the same way as the twist-2 . 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 can be repeated for as in the previous subsection. The factorized cross section reproduces the structure function up to the terms of order or .

At the level of individual contributions, the target subgraph introduces a non-zero 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 (the heavy-quark mass scheme). In particular, if the flavor-excitation coefficient function is modified by a term of order the twist-4 component of the structure function is modified by a term of order . The net change does not exceed the total error 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 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 . Therefore, the claim in Refs. Ball et al. (2016a, 2015a, b) 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 pairs in neutral-current DIS, and not as lone (anti)quarks, the rescaling in the S-ACOT- scheme Tung et al. (2002), requiring production in pairs only, approximates energy-momentum conservation better than its full ACOT counterpart that also tolerates production of single or 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.2.1, see also Aivazis et al. (1994b) and Guzzi et al. (2012).

The first line in Fig. 3 shows the lowest-order twist-4 contributions of order the remaining lines show some radiative contributions of order As before, a superscript in the parentheses indicates the order of the perturbative coefficient.

In either the 3- or 4-flavor scheme, we start by computing “flavor production” structure functions , such as or shown in Fig. 3, with standing for a or another double-parton initial state. Many more diagrams besides the ones shown arise at each order depending on the locations of the extra gluon attachments in the hard subgraph. The coefficient functions associated with twist-4 are derived by matching the perturbative coefficients order-by-order as in Eqs. (5) and (6).

For instance, at order , the double-convolution
integral over the gluon-pair
light-cone momentum fractions and scales as
where is at least
as large as or . In the limit ,
the contribution with a smaller hard scale
still survives. A part of it is resummed in the flavor-excitation
term , added across all
in order to obtain a smooth prediction for
.^{5}^{5}5The discontinuity of at the switching point
is very mild at NNLO and reduced with including higher
orders. Smoothness of is desirable for the convergence of
PDF fits.

The twist-4 remainder of that is not absorbed in may be of the same order as at relatively low . The remainder is given by where is found from the comparison of the coefficients in Eq. (6):

(7) |

The Feynman diagram for the “subtraction term” is shown in the first row of Fig. 3. It is obtained by inserting into the Feynman graph for 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 above the insertion by a simpler subgraph given by . Clearly, the remainder is process-dependent.

The next-order contribution with an added gluon line develops a logarithmic enhancement at ,

(8) |

which is resummed as a part of and . Again, the remainder that is not resummed must still be included in the full result, it takes the form of where

(9) |

stands for the infrared-safe part (with respect to light partons) of in the scheme Guzzi et al. (2012).

The rest of the coefficient functions can be computed along the same lines.

## Iii Models for the fitted charm

### iii.1 Overview

To recap the previous sections, a non-zero initial condition at for the “intrinsic charm PDF”, interpreted in the sense of the “fitted charm”, may be used to test for the power-suppressed charm scattering contribution of order , 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 for double-parton initial-states , as well as of twist-4 nonperturbative functions such as .

Various model estimates suggest a power-suppressed charm cross section of a modest size: of order of a fraction of the 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 Pumplin et al. (2007); Dulat et al. (2014). 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 at the initial scale will be considered. [ is set to be less than in all cases.] Besides the conventional CT14 model that sets , the other three models allow for of an arbitrary magnitude. In all models, the charm PDF is convoluted with the S-ACOT- coefficient functions with . It remains constant below the switching scale and is combined with the perturbative charm component at and evolved to 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.2.2: etc. The remainders are process-dependent and comparable to the convolutions at energies close to . 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 or , separation of the and terms presents an additional challenge. The experimental data in the CT14(HERA2) fits is selected with the cuts , so as to minimize sensitivity to the terms. This is usually sufficient to minimize the contributions below the PDF uncertainty from other sources. We examine the possibility of the impact of the terms on the best-fit in Sec. IV.5.

### iii.2 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 above
0.1 and satisfies with
for and with
for . The distributions for valence and quarks
fall into this broad category, as well as the “intrinsic” sea-quark
distributions that can be naturally generated in several ways
Pumplin (2006), e.g.,
for all flavors, nonperturbatively from a
Fock state in light-cone
Brodsky et al. (1980, 1981); Chang and Peng (2011); Blümlein (2016); Brodsky et al. (2015)
and meson-baryon models
Navarra et al. (1996); Paiva et al. (1998); Steffens et al. (1999); Hobbs et al. (2014); for and , from connected diagrams in lattice
QCD Liu et al. (2012).^{6}^{6}6In contrast to the light flavors,
in lattice QCD a charm PDF arises exclusively from
disconnected diagrams Liu (). This suggests
that and contributions in DIS are connected
to the hadron target by gluon insertions, in accord
with the physical picture of the
QCD factorization in Sec. II.1.

A sea-like component is usually monotonic in and satisfies for and for , with slightly above 1, and . This behavior is typical for the leading-power, or “extrinsic” production. For example, an (anti)quark PDF with this behavior originates from splittings in perturbative QCD, or from disconnected diagrams in lattice QCD (see Ref. Liu et al. (2012) for details). Even a missing next-to-next-to-next-to-leading order (N3LO) leading-power correction may produce a sea-like contribution at , 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 regions where it can be important (see the discussion of the EMC charm data in Sec. IV.6). 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 to be almost the same in size as or at and , 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 , while a negative-valued SEA form is not statistically distinguishable in the fit from a positive SEA form at a larger value. [The dependence of SEA fits on 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.

### iii.3 The charm distribution models in detail

We will now review these four models, whose distributions at GeV and 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 . The size of the preferred charm distribution at a given significantly depends on the values of the physical charm quark mass and QCD coupling strength . On the other hand, its dependence on the auxiliary theoretical scales of order , including the switching scale and the scale in the rescaling variable , cancels up to N3LO and thus is relatively weak; see a practical illustration in Fig. 1 of Gao et al. (2013). The net momentum fraction of the proton carried by charm quark starts off close to zero at and effectively saturates at high values at a level of approximately 2.5%, see Fig. 7.

ii) The approximate Brodsky-Hoyer-Peterson-Sakai (BHPS) model Brodsky et al. (1980, 1981) parametrizes the charm PDF at by a “valence-like” nonperturbative function

(10) |

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, 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 in Eq. 10. The parametrizations for and 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 . The parametrizations of this kind tend to have enhanced and ratios at , 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 is obtained by solving the BHPS model for the Fock state numerically and keeping the exact dependence on , and . This fit also includes small BHPS contributions to the and antiquarks generated from the and Fock states according to the same method. In the BHPS model, the quark distributions are determined by starting from a proton Fock state, where the probability differential for a quark to carry a momentum fraction is given by

(11) |

The standard BHPS result, used in ii), is given by letting and taking the limit to produce Eq. (10). However, Ref. Chang and Peng (2011) has shown that the solution that keeps the masses finite, including those of the light quarks, modifies the shape of , slightly shifting the peak to smaller . A similar conclusion was reached in Ref. Blümlein (2016), 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 .

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 and components, we obtain physically consistent ratios of the charm-quark and light-antiquark PDFs at large , cf. Fig. 9. We do not, however, include the BHPS contribution to the 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:

(12) |

This model is assumed with the SEA1 and SEA2 PDF sets from the two global fits distinguished by the value of normalization in Eq. 12.

Finally, the normalization coefficient in models ii)-iv) can be derived from the charm momentum fraction (first moment) at scale :

(13) |

By its definition, is evaluated at the initial scale . It is to be distinguished from the full charm momentum fraction at , which rapidly increases with because of the admixture of the twist-2 charm component.

## Iv Features of the CT14 intrinsic charm

### iv.1 Settings of the fits

The BHPS1, BHPS2, SEA1, and SEA2 parametrizations
are obtained by following the setup of the CT14 analysis
Dulat et al. (2016). BHPS3 is obtained with
the CT14HERA2 setup Hou et al. (2017a). 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 Arneodo et al. (1997)
was dropped in CT14HERA2, and the low- behavior of the strange
(anti)quarks was no longer tied to that of the and
antiquarks. This extra flexibility in of CT14HERA2
resulted in a reduction of over the entire
range relatively to CT14. This feature has
potential implications for the models of with a sea-like
behavior. In some exploratory fits, we include the EMC data Aubert et al. (1983)
on semiinclusive DIS charm production, while in the other fits we
examine sensitivity on the input pole charm
mass.^{7}^{7}7CTEQ-TEA fits can also take a charm mass, rather than
the pole mass as the input Gao et al. (2013), with similar
conclusions.

The PDFs for light partons are parametrized at an initial scale slightly below GeV, with the exception of the study of the dependence, in which it was more convenient to start at a lower initial scale GeV. For all models, the QCD coupling constant is set to , compatible with the world average value Patrignani et al. (2016) , as in the standard CT PDF fits. The PDFs are evolved at NNLO with the HOPPET code Salam and Rojo (2009). NLO ApplGrid Carli et al. (2010) and FastNLO Wobisch et al. (2011) interpolation interfaces, combined with NNLO/NLO factor look-up tables, were utilized for fast estimation of some NNLO cross sections.

### iv.2 Dependence on the charm momentum fraction

In the models in Sec. III.3, the magnitude of is controlled by its normalization , correlated uniquely with the net momentum fraction of defined in Eq. (13). The choice of the 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 , we examine the goodness-of-fit function

(14) |

constructed in the CT14 method from the global and a “tier-2” statistical penalty Dulat et al. (2016). It is convenient to compare each fit with an to the “null-hypothesis” fit obtained assuming . Thus, we start by computing

(15) |

where and are given for and , respectively, at 50 values of and default GeV. We plot the resulting 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 without the tier-2 penalty , and a solid one for that includes according to Eq. (14). The function estimates the global quality of fit and is equal to the sum of contributions from all experiments and theoretical constraints. A non-negative “Tier-2” penalty is added to to quantify agreement with each individual experiment Dulat et al. (2014); Gao et al. (2014). Being negligible in good fits, grows very rapidly when some experiment turns out to be inconsistent with theory. The net effect of is to quickly increase the full if an inconsistency with some experiment occurs, even when 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 , however, is observed for the BHPS fits, roughly at , both in the CT14 and CT14HERA2 frameworks.

The significance of this reduction and the upper limit on depends on the assumed criterion. In CTEQ practice, a set of PDFs with smaller (larger) than 100 units is deemed to be accepted (disfavored) at about 90% C.L. Thus, a reduction of 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 in the CT14 and CT14HERA2 analyses at the 90% C.L.:

(16) |

In keeping with the previous analysis of Ref. Dulat et al. (2014), we define specific fits with particular choices of for both examined models. The fits BHPS1 and SEA1 correspond to , while BHPS2 has and SEA2 has . 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 , 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 and quarks.

The results in Fig. 5 are compatible with
the findings of the previous CT10 NNLO IC analysis Dulat et al. (2014).
In particular, comparing to CT14 in the left frame of
Fig. 5 and to Fig. 2 in
Ref. Dulat et al. (2014),^{8}^{8}8In
Ref. Dulat et al. (2014), , , and are
denoted by , , and . we see
that the minimum in 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 of the SEA model rises more rapidly with increasing than it does in the comparable CT14 fit. This is due to the greater flexibility in the low- behavior of the strange-quark distribution in the CT14HERA2 framework discussed previously. More freedom reduces at low and thus increases and at the same . In the CT14 fit with the SEA charm component, the minimum is at , and it is largely washed out in the CT14HERA2 case. The for SEA grows faster for CT14HERA2 compared to CT14: at it is higher by about 40 units in Fig. 5(right) relatively to Fig. 5(left).

The reduction in for the NNLO BHPS fits at = 0.01, relatively to the fit with , thus remains a persistent feature of the CT10, CT14, and CT14HERA2 analyses. While the 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.1? It will be discussed in Sec. IV.5 that is reduced primarily in a few fixed-target experiments (the 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.

### iv.3 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 is very sensitive to the assumed charm mass.

Distinct from the auxiliary QCD mass parameters – , , 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 , not to the auxiliary parameters of order Gao et al. (2013). The dependence remains mild, the values in the range GeV are broadly consistent with the CT14 data.

Exploratory fits testing the dependence of on , for a selection of pole masses, 1.1, 1.2, 1.3, 1.4, 1.5 GeV, are illustrated by Fig. 6. The general setup of these scans follows the fits to the CT14 data. To access the masses below 1.3 GeV, we reduced the initial scale to 1 GeV and examined alternative forms for the gluon PDF parametrization, because DIS charm production is sensitive to the gluon PDF . Dependence of for CT14 NNLO on for two representative gluon parametrizations at 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 GeV, is constrained to be positive at all ; while for “gluon 2”, it is allowed to be negative at the smallest and , 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 measured by the H1 Collaboration Aaron et al. (2011). The more flexible “gluon 2” parametrization results in a marginally better with respect to the nominal CT14, or “gluon 1”, at a slightly lower GeV, and with a large uncertainty. This best-fit value in this range is consistent with the previously observed tendency of the DIS data to prefer smaller masses at , e.g., GeV obtained in the CT10 setup Gao et al. (2013).

The two lower insets of Fig. 6 illustrate the variations in , with the more flexible “gluon 2” parametrization, when the IC component is included for five values of . The circles on the curves mark the minima; the thin vertical lines indicate the exclusion limits on for each value.

For the BHPS model in the left inset, the position of the minimum is relatively stable as is varied, while the upper limit on decreases to 1.9% as increases. The overall conclusion is that the preferred at scale is not strongly sensitive to the variations of in the case of the BHPS parametrizations. On the other hand, as we will see in a moment, the total momentum fraction at scales above is sensitive to due to the growing perturbative charm component.

The situation is very different for the SEA model shown in Fig. 6 (right), where the dependence on is more pronounced. In this case, develops a pronounced minimum for GeV, while the minimum totally disappears, and is totally excluded, for GeV.

This can be understood as follows: when increases, the twist-2 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 , which enhances the scattering contribution. An opposite effect occurs when decreases (i.e., less phase-space suppression for fusion, a smaller intrinsic charm momentum allowance in scattering). But the and