New CTEQ global analysis of quantum chromodynamics with high-precision data from the LHC

We present the new parton distribution functions (PDFs) from the CTEQ-TEA collaboration, obtained using a wide variety of high-precision Large Hadron Collider (LHC) data, in addition to the combined HERA I+II deep-inelastic scattering data set, along with the data sets present in the CT14 global QCD analysis. New LHC measurements in single-inclusive jet production with the full rapidity coverage, as well as production of Drell-Yan pairs, top-quark pairs, and high-$p_T$ $Z$ bosons, are included to achieve the greatest sensitivity to the PDFs. The parton distributions are determined at NLO and NNLO, with each of these PDFs accompanied by error sets determined using the Hessian method. Fast PDF survey techniques, based on the Hessian representation and the Lagrange Multiplier method, are used to quantify the preference of each data set to quantities such as $\alpha_s(m_Z)$, and the gluon and strange quark distributions. We designate the main resulting PDF set as CT18. The ATLAS 7 TeV precision $W/Z$ data are not included in CT18, due to their tension with other data sets in the global fit. Alternate PDF sets are generated including the ATLAS precision 7 TeV $W/Z$ data (CT18A), a new scale choice for low-$x$ DIS data together with a slightly higher choice for the charm mass (CT18X), or all of the above (CT18Z). Theoretical calculations of standard candle cross sections at the LHC (such as the $gg$ fusion Higgs boson cross section) are presented.


INTRODUCTION
With an accumulated data sample of over 140 fb −1 at the 13 TeV run for both ATLAS and CMS collaborations, the Large Hadron Collider (LHC) has entered an era of precision physics. The experimental precision has been matched by improvements to the theoretical predictions, with a number of collider processes now available at the next-to-next-to-leading order (NNLO) in the QCD coupling strength. Such precision is necessary for rigorous tests of the Standard Model (SM) and in searches for signs of physics beyond the Standard Model (BSM), as there have been no 'smoking-gun' signs of BSM physics to date. Precise predictions in QCD theory require correspondingly precise parton distribution functions (PDFs) [1][2][3][4][5][6][7][8], which in turn warrant advances in interpreting LHC experiments to extract important information about the SM, and possibly, BSM physics.
To this end, we present a new family of CTEQ-TEA parton distribution functions, designated as CT18. These PDFs are produced at both next-to-leading order (NLO) and NNLO in the QCD coupling constant, α s . The CT18 PDFs update those of the CT14 family presented in Ref. [1]. In the new analysis, we include a variety of new LHC data, at the center-of-mass energies of 7 and 8 TeV, on production of single-inclusive jets, W/Z bosons, and top-antitop quark pairs, obtained by the ATLAS, CMS and LHCb collaborations. At the same time, the update retains crucial "legacy" data from the previous CT global QCD analyses, such as the HERA I+II combined data on deep-inelastic scattering (DIS) and measurements in fixed-target experiments and at the Fermilab Tevatron pp collider. Measurements of processes in similar kinematic regions, by both ATLAS and CMS, allow crucial cross-checks of the data. Measurements by LHCb often allow extrapolations into new kinematic regions not covered by the other experiments. Some processes, such as tt production, allow for the measurement of multiple observables that provide similar information for the determination of PDFs. In addition to the PDFs themselves, we also present relevant PDF luminosities, and predictions with uncertainties for standard candle cross sections at the LHC.
The goal of the CT18 analysis is to include as wide a kinematic range for each measurement as possible, while still achieving reasonable agreement between data and theory. For the ATLAS 7 TeV jet data [9], for example, all rapidity intervals cannot be simultaneously used without the introduction of systematic error decorrelations provided by the ATLAS 6 collaboration [10]. Even with that decorrelation, the resultant χ 2 for the new jet experiments is not optimal, resulting in less effective PDF constraints. Inclusive cross section measurements for jet production have been carried out for two different jet-radius values, R, by both ATLAS and CMS. For both experiments, we have chosen the data with the larger R-value, for which the NNLO (fixed order) prediction should have a higher accuracy.
We evaluate the jet cross section predictions using a QCD scale of inclusive jet transverse momentum Q = p jet T , consistent with past usage at NLO. The result is largely consistent with similar evaluations using Q = H T [11][12][13].
Theoretical predictions for comparison to the data used in the global fit have been carried out at NNLO, either indirectly through the use of fast interpolation tables such as fastNLO [14,15] and ApplGrid [16], together with NNLO/NLO K-factors, or directly (for top-quark related observables) through the use of fastNNLO grids [17,18].
In an ideal world all such data sets would perfectly be compatible with each other, but differences are observed that do result in some tension between data sets and pulls in opposite directions. One of the crucial aspects of carrying out a global PDF analysis is dealing with data sets that add some tension to the fits, while preserving the ability of the combined data set to improve on the existing constraints on the PDFs. In some cases, a data set may be in such tension as to require its removal from the global analysis, or its inclusion only in a separate iteration of the new PDF set.
In this paper, we will describe how the high-precision ATLAS 7 TeV W/Z rapidity distributions, which, as we find, favor an increase of the strange quark distribution at low x, require such special treatment. In particular, while other PDF-analysis groups (e.g., MMHT, see Ref. [19]) have noted that these data can be fitted with χ 2 /d.o.f. that is comparable to the CT18 one, we find that such χ 2 reflects systematic tensions with many of the other data in our global analysis if this information is included in CT18. Furthermore, the standard Hessian profiling technique used by the experimental collaborations significantly underestimates the minimal χ 2 that can be reached for the ATLAS 7 TeV W/Z data when they are included in the CT18 fit. We therefore treat these measurements separately in an alternative fit, CT18A, introduced in Sec. II C. In another variant, CT18X, a special scale is used for the calculation of low-x DIS cross sections; the scale mimics the impact of low-x resummation.
Both modifications cause an increase in the low-x quark and gluon distributions. Finally, these two fit variants are amalgamated into a combined alternative fit, CT18Z, canvassed in Appendix A, and to which we compare results throughout this article.
Our current global analyses are carried out in four stages. First, PDFSense [20,21], a program for a rapid survey of QCD data using the Hessian approach [22,23], is used to select the data sets that are expected to have the greatest impact on the global PDF sets.
This selection takes into account both the correlation between the data and specific PDFs, in a given x range, and the sensitivity of the data to those PDFs. The sensitivity depends on the magnitudes of the statistical and systematic errors (and the correlations of the latter).
For example, both the collider inclusive jet data and the top-quark data have a strong correlation with the high-x gluon, but the inclusive jet data has a larger sensitivity due to a much larger number of data points. Next, ePump [24,25] is used to quickly examine the quantitative impact of each selected data set, within the Hessian approximation. Third, the full global PDF fit is carried out using all such data sets. Recent enhancements to the CT global analysis code have greatly improved the speed of the calculations. Lastly, the impact of key data sets on certain PDFs at specific kinematic points of interest, as well as on the value of α s (M Z ), is assessed using the Lagrange Multiplier (LM) method [26].
In order to minimize any parametrization bias, we have tested different parametrizations for CT18: e.g., using a more flexible parametrization for the strange quark PDF. In some kinematic regions, there are fewer constraints from the data on certain PDFs. In these cases, LM constraints have been applied to limit those PDFs to physically reasonable values.
Our paper is organized as follows. Section II begins with an executive summary of the key stages and results of the CT18 global analysis. It continues with an overview of the chosen experimental data and alternative fits (CT18Z, CT18A, and CT18X) in the CT18 release.
In Sec. III we detail theoretical/computational updates to the CT fitting methodology and details for specific process-dependent calculations. Sec. IV presents the main results obtained in CT18 -the fitted PDFs as functions of x and Q, determinations of QCD parameters (α s , m c ), calculated parton luminosities, and various PDF moments and sum rules. This is the heart of the paper, it will be of interest to a broad group of researchers who will use the PDFs for theoretical predictions at LHC experiments. Sec. V describes the ability of CT18 to provide a successful theoretical description of the fitted data. In addition to characterizing the fit of individual data sets, in Sec. VI we also compute the various standard candle quantities of relevance to LHC phenomenology, for instance, Higgs boson production cross sections at 13 and 14 TeV, and various correlations among electroweak boson and top-8 quark pair production cross sections. In Sec. VII, we discuss the broader implications of this work and highlight our main conclusions. Several appendices present a number of important supporting details. In Appendix A, we review the CT18Z and other alternative fits, including descriptions of various data sets admitted into these separate analyses. A number of more formal details related to our likelihood functions and relations among covariance matrices are summarized in Appendix B. Appendix C presents the analytical fitting form adopted in CT18 and values of fitted parameters. Appendix D presents a number of technical advances in the CT fitting framework, including code parallelization, while Appendix E enumerates the decorrelation models utilized in fitting the newly included inclusive jet data from the LHC. Lastly, in Appendix F, we present the results of a short study based on Hessian profiling methods to assess the impact of the 7 TeV W/Z production data taken by ATLAS.

II. OVERVIEW OF THE CT18 GLOBAL QCD ANALYSIS
A. Executive Summary

Input experimental data and final PDF ensembles
The CT18 analysis updates the widely used CT14 PDF sets [1] by applying NNLO and NLO global fits to an expanded set of experimental measurements that include high-statistics data from the ep collider HERA and the LHC. The CT18 experimental data set includes high-statistics measurements from ATLAS, CMS, and LHCb on production of inclusive jets, W/Z bosons, and top-quark pairs, while it retains the crucial legacy data, such as measurements from the Tevatron and the HERA Run I and Run II combined data. The kinematic distribution of the data points used in CT18 is shown in Fig. 1 in a space of parton momentum fraction, x, and QCD factorization scale, denoted here as Q. As has been true for global PDF fits for some time, the data included cover a large kinematic range, both in x and Q. We review implementation of the new data sets in the CT18 global fit, the associated physics issues that affect the resulting PDFs, and a wide class of QCD predictions based on them.
By 2018, the LHC collaborations published about three dozen experimental data sets that can potentially constrain the CT PDFs. In light of the unprecedented precision reached in some measurements, the latest LHC data must be analyzed using next-to-next-to-leading     [20]. Also shown are the ATLAS 7 TeV W/Z production data (ID=248), labeled ATL7WZ'12, fitted in CT18Z. order (NNLO) theoretical predictions in perturbative QCD. The fitted PDFs we obtain in this analysis are plotted in Fig. 2, which displays in the upper panels the CT18 PDFs at two widely-separated scales, Q = 2 and 100 GeV (on the left and right, respectively). In the lower panels, we show the corresponding PDFs found in our amalgamated alternative analysis, CT18Z.
Among the four ensembles of PDFs, CT18, A, X, and Z, the CT18 and CT18Z ensembles are the most dissimilar in terms of the shapes of PDFs in the LHC kinematic range.
A number of features are apparent by comparing the CT18 and CT18Z families of fitted PDFs, including changes in the x dependence of the fitted strangeness distribution s(x, Q) and associated effects in the extracted uncertainty for CT18Z. At the same time, for CT18, we obtain modest improvements in the precision for the gluon density g(x, Q), as compared GeV for u, u, d, d, s = s, and g. Lower panels: The analogous curves, but obtained for CT18Z. In all instances, the gluon PDF has been scaled down as g(x, Q)/5. The charm distribution, c(x, Q), which is perturbatively generated by evolving from Q 0 = 1.3 and 1.4 GeV, respectively, in CT18 and CT18Z, is also shown.
to CT14, following the inclusion of the LHC Run-1 data discussed below. For CT18Z, however, we obtain a somewhat enlarged uncertainty for the gluon and perturbatively-generated charm PDFs, especially at the lowest values of x < 10 −3 , due to modified treatment of the DIS data as described in Sec. II C and App. A. These final PDFs depend on numerous systematic factors in the experimental data. Scrupulous examination of the systematic effects was essential for trustworthy estimates of PDF uncertainties, and the scope of numerical computations also needed to be expanded.  Even in the LHC era, DIS data from the ep collider HERA provide the dominant constraints on the CT18 PDFs. This dominance can be established using the ePump and PDFSense statistical techniques reviewed below. CT18 implements the final ("combined") data set from DIS at HERA Run-I and Run-II [27] that supersedes the HERA Run-I only data set [28] used in CT14 [1]. A transitional PDF set, CT14 HERAII , was released based on fitting the final HERA data [29]. We found fair overall agreement of the HERA I+II data with both CT14 and CT14 HERAII PDFs, and that both PDF ensembles describe equally well the non-HERA data included in our global analysis. At the same time, we observed some disagreement ("statistical tension") between the e + p and e − p DIS cross sections of the HERA I+II data set. We determined that, at the moment, no plausible explanation could be provided to describe the full pattern of these tensions, as they are distributed across the whole accessible range of Bjorken x and lepton-proton momentum transfer Q at HERA.
Extending these studies using the CT18 fit, we have investigated the impact of the choice of QCD scales on inclusive DIS data in the small-x region, as will be explained later in Sec. II C.
We find that the quality of fit to HERA data is improved by about 50 units by evaluating the NNLO theoretical cross sections in DIS with a special x-dependent factorization scale, µ F,x , introduced in Section II C. Fig. 3 (left) shows the changes in the candidate CT18 PDFs obtained by fitting the DIS data sets with the factorization scale µ F,x , as compared to the CT18 PDFs with the nominal scale µ F = Q. With the scale µ F,x , we observe reduced u and d (anti-)quark PDFs and increased gluon and strangeness PDFs at x < 10 −2 , as compared to the nominal CT18 fit, with some compensating changes occurring in the same PDFs in the unconstrained region x > 0.5 in order to satisfy the valence and momentum sum rules.
The right panel of Fig. 3 shows the χ 2 /N pt values (divided by the number N pt of experimental data points) for four HERA data sets (inclusive neutral and charged current DIS [27], reduced charm, bottom production cross sections, and H1 longitudinal function [30]) in the fits as a function of the statistical weight w of the HERA I+II inclusive DIS data set [27]. The default CT18Z fit correspond to w = 1; with w = 10, the CT18Z fit increasingly behaves as a HERA-only fit. We see that, with the scale µ 2 F,x and w = 10, χ 2 /N pt for the inclusive DIS data set improves almost to the levels observed in the "resummed" HERA-only fits without intrinsic charm [31,32]. The quality of the fit to the charm semi-inclusive DIS (SIDIS) cross section and H1 F L also improves.

Selection of new LHC experiments
When selecting the most promising LHC experiments for the CT18 fit, we had to address a recurrent challenge -the presence of statistical tensions among various (sub)sets of the latest experimental data from HERA, LHC, and the Tevatron. The quickly improving precision of the collider data reveals previously irrelevant anomalies either in the experiment or theory. These anomalies are revealed by applying strong goodness-of-fit tests [33]. Figure 4 illustrates the degree of tensions using a representation based on the effective Gaussian vari- [34] constructed from the χ 2 values and numbers of data points N E for individual data sets E. In a high-quality fit, the probability distribution for S E must be approximately a standard normal distribution (with a unit half-width). In the global fits by CTEQ-TEA and external groups, we rather observe wider S E distributions as in Fig. 4

Advances in fitting methodology
The CTEQ-TEA group resolved these challenges through a multi-pronged effort. We developed two programs for fast preliminary analysis to identify the eligible experimental data sets for the global fit. The PDFSense program [20] was developed at Southern Methodist University (SMU) to predict quantitatively, and before doing the fit, which data sets will have an impact on the global PDF fit. The ePump program [24]  The four families of general-purpose CT18 PDFs presented in this paper are most suitable for computations of central predictions and estimations of PDF uncertainties in accord with the guidelines on the usage of PDFs published by the PDF4LHC group [52]. The PDF uncertainties are constructed at the 90% probability level based on two tiers of criteria as in the CT14 global analysis [1]. These PDFs are obtained by assuming a world-average QCD coupling constant, α s (M Z ) = 0.118 [53]. The combined PDF+α s uncertainty can be computed using the special α s series of the PDFs for each family by adding the PDF and α s uncertainties in quadrature, as explained in Ref. [54].
As we have done with the previous families of the CTEQ-TEA PDFs, we intend to produce specialized PDF sets to complement the current general-purpose PDF ensembles.
In particular, CT18 updates on the PDFs with an asymmetric-strangeness PDF [55], photon PDF [56] and fitted charm [57] will be provided in future publications.

B. Experimental data sets fitted in CT18
The CT18 global analysis starts with the data set baseline of CT14 HERAII [29] and adds the latest precision LHC results. The experiments in the CT14 HERAII baseline are listed in Table I, while the new LHC data sets included in the CT18(Z) fit are shown in Table II.
Tables I and II also include information on the number of data points, χ 2 , and the effective Gaussian variable, S E , for each individual data set appearing in the global fit. Most of the data sets are included in all four PDF ensembles; we will identify differences between the specific selections as they arise.

Charting sensitivity of new data sets to the PDFs
As discussed in Secs. II A 3 and II A 4, we employed new methods, based on the Hessian correlation between "data minus theory" residuals [22,58,59], and on a more informative "sensitivity variable" described in Ref. [20], to determine quantitatively a hierarchy of impact of data on the global fit, and on specific cross sections.
In particular, we used the PDFSense framework to determine which new data sets would be most important for the CT18 update. As a demonstration of this, we use it here to predict which data sets have the most impact on one of the most crucial predictions at the LHC, that of the Higgs boson cross section (σ H ) through the gg fusion process (at √ s = 14 TeV). Often, to get an indication which data sets will have the most impact on such cross sections, one examines the Pearson correlations between the experimental data points and the gluon distribution in the kinematic region responsible for Higgs boson production. The left-hand side of Fig. 5 shows the data points with the highest absolute correlations |C f | (defined by the statistical residuals as in Ref. [20]) with the Higgs boson cross section at 14 TeV.
By this measure, there may be a number of high-impact data sets, notably HERA neutral current DIS, LHC and Tevatron jet production, and HERA charm quark production. The correlations, however, do not reflect the experimental uncertainties of the data points: an experimental cross section could be highly correlated with the gluon distribution in the x range responsible for Higgs boson production and still not provide much of a constraint on the Higgs boson cross section if the experimental uncertainties are large. Conversely, an   [20]. Right: A similar assessment of the CT18 candidate data, but computed on the basis of the sensitivity, |S f |(x i , Q i ). In both panels, a highlighting cut has been imposed to draw attention to the ∼ 300 highest-impact points according to each metric. experimental cross section that might not have as large a correlation, but which has smaller (statistical and systematic) uncertainties, may provide a stronger constraint.
The level of constraint is thus better predicted by the sensitivity variable S f , defined in Ref. [20]. The experimental data points used in the CT18 global fit that have the highest absolute sensitivity |S f | to the PDF dependence of σ H at 14 TeV are shown on the right-hand side of Fig. 5. More data points (from a larger number of experiments) have high sensitivity than those identified by high correlation. In addition to the DIS data from HERA I+II, there are also contributions from the fixed-target DIS experiments, as well as measurements from the LHC.
As will be shown in Sec. V A 1, using Lagrange Multiplier scans, the HERA I+II data still dominate the determination of the gluon distribution in the range sensitive to Higgs boson production at the LHC, because of its small statistical and systematic errors. Because of the continuing influence of these older data sets, we will find that the reduction of the PDF uncertainty for Higgs boson production is less significant in CT18 than in CT14. In addition, tension between some of the most sensitive data sets limits the reduction on the uncertainty of the Higgs cross section. These effects are explored in detail in Sec. V A 1.
In order to corroborate the findings of PDFSense, we added the new LHC data sets one-byone on top of the CT14 HERAII NNLO baseline data using the software package ePump [24,25], for updating and optimizing error PDFs in the Hessian approach. The ePump package can quickly assess the impact of each new data set on the PDFs without performing a full global fit. Furthermore, Lagrange multiplier studies have been carried out to examine constraints of specific data sets on PDF distributions and to explore the impact of PDF parametrizations at large or small x, and (in some cases) to highlight tension between experiments. Finally, a parallelization of the global PDF fitting package has been implemented in order to allow for faster turn-around time.
The final CT18(Z) data ensemble contains a total of N pt = 3681(3493) data points and results in χ 2 /N pt = 1.17 (1.19) at NNLO. We now discuss some of the new data sets included in the CT18 analysis, and highlight the differences in the CT18Z fit.

Baseline data sets
The CT18 global analysis inherits from CT14HERA2 a number of precision non-LHC experiments listed in Table 1. Among those, the HERA I+II DIS data set provides the most significant constraints, followed by a group of fixed-target neutral-current DIS experiments: BCDMS, NMC, and CCFR. Similarly, a number of neutrino DIS measurements have previously been included and provide valuable constraints on sea (anti-)quarks. Among them, we find that the single-nucleon structure functions F p 2 and F p 3 extracted from CDHSW data on neutrino-iron deep inelastic scattering exhibit a preference for a harder gluon PDF at x 0.1, compared to CCFR and other experiments, cf. Fig. 20. This well-known behavior reflects larger logarithmic slopes of F p 2 and xF p 3 measured by CDHSW, as compared to the analogous CCFR measurements [60], which in turn may reflect differences in the energy calibration and resolution smearing between the two experiments [61]. Thus, to help obtain a softer large-x gluon behavior, as being favored by recent LHC data, we exclude the CDHSW F 2 and xF 3 data sets from the CT18Z analysis, while including these sets in the rest of the CT18 PDF ensembles.
We continue to include a variety of lepton pair production measurements from the Teva-tron and fixed-target experiments, as summarized in Table I. The low-statistics data on W/Z production at LHCb 7 TeV [62] and ATLAS, CMS 7 TeV jet production [63,64] are replaced in the CT18 analyses by more recent measurements, as summarized in the next section.

LHC precision data from W/Z vector boson production
The CT18(Z) global analysis uses W/Z vector boson production data from LHC Run-I, including measurements from the ATLAS, CMS and LHCb collaborations.
The measurements from ATLAS included in the fit are: • The √ s = 7 TeV W/Z combined cross section measurements [35] (ID=248) with 4.6 fb −1 of integrated luminosity. The ATLAS group has performed 7 measurements with a total of 61 data points: distributions in the pseudorapidity of charged lepton in W + (11 points) and W − (11 points) production; rapidity of lepton pairs for low-mass Drell-Yan (DY) process in the central region (6 points); Z-peak DY process in the central (12 points) and forward (9 points) regions; high-mass DY process in the central (6 points) and forward (6 points) regions. In the published fits, we include 3 measurements: W + , W − and Z-peak central DY production (34 points in total). These data are used only to fit the CT18A and CT18Z PDFs, but not the CT18 and CT18X PDFs. Other data are ignored due to the sizable EW corrections and/or photon-induced contribution (γγ → l + l − ), as discussed in Sec. V C.
• The √ s = 8 TeV distribution of transverse momentum p T of lepton pairs in the Z/γ * production (ID=253) [95] with 20.3 fb −1 of integrated luminosity. The ATLAS collaboration measured the p ll T distribution up to 900 GeV for the lepton pairs in the invariant mass range 12 < M ll < 150 GeV. Meanwhile, the experimentalists presented both the normalized and absolute cross sections for the singly differential distribution dσ/dp ll T and doubly differential distribution d 2 σ/(dp ll T dy ll ). To select the cleanest and most sensitive data for the CT18 fits, we only include 3 invariant-mass bins around   data at M ll < 46 GeV, for which the kinematic cut p l T > 20 GeV restricts the cross section to be coming predominantly from the region p ll T M ll , where the higher-order corrections beyond the current O(α 3 s ) calculation are significant. We fit neither the normalized p ll T distributions, as they introduce artificial interdependence between the particle rates in the disparate p ll T regions through their shared overall normalization, nor the doubly differential distributions, which are not very sensitive. In total, we include 27 data points in the pair's transverse momentum region 45 ≤ p ll T ≤ 150 GeV, where the fixed-order NNLO cross section is most reliable. The data at lower p ll T and higher p ll T regions are excluded because of contributions from small-p T resummation and electroweak corrections, as discussed in Sec. V C.
For CMS, measurements of the charge asymmetry for inclusive W ± production at √ s = 8 TeV [93] (ID=249) are included, with 18.8 fb −1 of integrated luminosity. These consist of 11 bins of muon pseudo-rapidity (over the range 0 ≤ |η µ | ≤ 2.4) with p µ T ≥ 25 GeV. The correlated systematic errors are implemented using a decomposition of the covariance matrix to convert it to the correlation matrix representation according to the procedure described in Appendix B. The same decomposition method will also be applied to the LHCb W/Z experiments (Exp ID 245 and 250) to convert the published covariance matrix to a correlation matrix. We have explicitly verified their equivalence using ePump, which is capable of operating with either covariance matrix or correlation matrix [25] representation.
In CT18, we include three experimental data sets published by LHCb: • The √ s = 7 TeV W/Z forward rapidity cross section measurements [91] (ID=245), • The √ s = 8 TeV Z → e + e − cross section measurements [92] (ID=246) at forward rapidity, with 2.0 fb −1 of integrated luminosity, consist of 17 bins of Z-boson rapidity (2.0 ≤ y Z ≤ 4.25). The luminosity uncertainty is taken to be fully correlated, but the other correlated uncertainties are simply added in quadrature. We have used ePump to confirm that this approximation yields similar updated PDFs as those obtained from using the covariance matrix representation.
As in the √ s = 7 TeV case, the correlated systematic errors are included by converting the covariance matrix into the correlation matrix representation. The beam energy and luminosity uncertainties are taken to be fully correlated between the cross-section measurements.

LHC inclusive jet production
For CMS, double-differential cross section measurements, d 2 σ/(dp T dy), for jet production at both √ s = 7 TeV and 8 TeV, are used. We use the larger of the two jet radii (R = 0.7) [96].

23
The data sets consist of 5 fb −1 of integrated luminosity at In addition to the systematic error information provided in the HEPData files, jet energy corrections (JEC) in the CMS 7 TeV data have been decorrelated according to the procedure in Ref. [100]. In particular the JEC2 ("e05") and an additional CMS-advocated decorrelation for |y| > 2.5, have been implemented [101]. These decorrelations improve the ability to fit the data.
For ATLAS, we again use the larger of the two jet radii (R = 0.6). Inclusive jet cross section measurements at √ s = 7 TeV with R = 0.6 and 4.5 fb −1 of integrated luminosity [9] (Exp. ID=544) are included in the global fit. This data set contains six bins covering the rapidity range 0.0 ≤ y ≤ 3.0, with a total of 140 data points in the 74 ≤ p T ≤ 1992 GeV range. This data set replaces the previous data set [63] which contains 37 pb −1 data.
Following the prescription given in Ref. [10], two jet energy scale (JES) uncertainties have been decorrelated in the ATLAS 7 TeV jet data, namely: MJB (fragmentation) ("jes16"), and flavor response ("jes62"). The decorrelation procedure reduces the χ 2 value by approximately 92 units. The total contribution to the χ 2 from the systematic error shifts is 28 (for 74 correlated systematic errors). Only one of the systematic error sources requires a shift greater than 2 standard deviations. A further improvement of 52 units is obtained by including a 0.5% theory error to account for statistical noise associated with the Monte Carlo calculations of the needed NNLO/NLO K-factors [11][12][13] in the NNLO fit. More details on the treatment of the ATLAS inclusive jet data are provided in Appendix E.

LHC top-quark pair production
ATLAS and CMS have measured top-quark pair production differential cross sections as a function of the top-(anti)quark transverse momentum p T,t , invariant mass m tt , rapidity of the pair y tt , transverse momentum of the top-quark pair p T,tt , and top-quark rapidity y t , individually for ATLAS, and in pairs for CMS. The individual impacts of the single differential tt cross section measurements have been analyzed, first by using the PDFSense sensitivity framework of Ref. [20], and second in separate fits via ePump in Refs. [25,102].
There is some tension between the observables, with respect to the gluon distribution that each prefers. We have made a selection of the top-quark pair production measurements to be included in the CT18 global analysis based on their compatibility within the fit.
Difficulties in fitting simultaneously p T,t , m tt , y t , and y tt distributions at 8 TeV were also found in [25,102,103].
By making use of the statistical correlations between observables in ATLAS, it is possible to include more than one of the observables in the global fit. We have chosen the absolute differential cross sections of the top-p T spectrum dσ/dp T,t , and the invariant mass dσ/dm tt , at √ s = 8 TeV with 20.3 fb −1 of integrated luminosity (Exp. ID=580) [99], based on the recommendation from ATLAS. 1 For CMS, we have chosen the normalized double differential The two ATLAS measurements are combined into one single data set which includes the full phase-space absolute differential cross-sections after the combination of the e+jets and µ+jets channels for the p T,t and m tt distributions with statistical correlations. Both of these distributions are fitted together by decorrelating one of the systematic uncertainties relative to the parton shower (PS) [104]. The QCD theoretical predictions at NNLO for these observables are obtained by using fastNNLO tables provided in Refs. [17,18].
The observed effect of these data sets on the CT18 PDFs is modest if the tt data are included together with the Tevatron and LHC jet production. The impact on the gluon PDF is consistent with the jet data, but the jet data provide stronger constraints, consistent with the larger number of data points, the wider kinematic range, and the relatively small statistical and systematic errors of the jet data.
In the course of the CT18 analysis, CMS measurements of top-quark pair production differential cross sections at 13 TeV were published [105], and bin-by-bin data correlations were made available on the HEPData repository. While these measurements are not currently included in the CT18 global fit, their agreement with theory is discussed later in Sec. VI.

Other LHC measurements not included in the CT18 fits
Besides the CT18(Z) data ensemble, we have carefully investigated the impact on the CT global fit of several other high-luminosity measurements from LHC Run-I. In certain cases we observed either no significant impact or substantial tensions with the CT18(Z) baseline.
The following vector boson production data were examined using PDFSense, ePump, or full fits, but not included in the final CT18(Z) global analysis: • Difficulties were encountered in obtaining a good agreement between theory and the ATLAS √ s = 7 TeV Z-boson transverse momentum distribution (p ll T ) data with 4.7 fb −1 of integrated luminosity [106]. The subset of these data with p ll T ∼ m ll (in the kinematic region most amenable to a fixed-order calculation) has rendered unacceptably high χ 2 values for various combinations of the renormalization and factorization scales that we have tried. No significant constraints could be ascribed to the CMS √ s = 8 TeV p ll T and y ll distributions with 19.7 fb −1 of integrated luminosity [107] in the Z peak kinematic region, and to the CMS √ s = 8 TeV normalized W p T and Z p T spectra with 18.4 pb −1 [108]. In the case of comparing the CMS double-differential distributions in (p T , y), we observed a large discrepancy between theory and data in the last rapidity bin. Since both ATLAS and CMS measurements are presented as normalized p ll T distributions of lepton pairs, it was not clear how to consistently compare to data in a limited range p ll T ∼ m ll when the normalization of data points was dependent on the cross section outside of the fitted range.
• No substantial changes in the candidate PDFs were observed after including either the single-or double-differential distributions, dσ/dQ or d 2 σ/(dQ dy), of the ATLAS √ s = 8 TeV Drell-Yan cross section measurements at 116 ≤ Q ≤ 1500 GeV and 0 ≤ y Z ≤ 2.5 with 20.3 fb −1 of integrated luminosity [109]. These high-mass data are impacted by non-negligible EW corrections and photon-induced (PI) dilepton production, the point that is further addressed in Sec. V C. For the same reason, we do not include the data of ATLAS 7 TeV high-mass Drell-Yan production with 4.7 fb −1 of integrated luminosity [110].
• The low-mass Drell-Yan data by the ATLAS collaboration at 7 TeV [111] were also explored and found to have no significant impact on the PDFs.   In light of the above discussion, in addition to the primary fit CT18, we present a number of alternative QCD global fits in order to disentangle the subtleties involved in describing the HERA and LHC Run-1 data. These additional fits were explored in parallel with CT18, by making alternative choices for data selection and theoretical calculations. The key differences among these fits are listed in Table III. Their predictions will be compared in Appendix A.
(i) CT18X differs from CT18 in adopting an alternate scale choice for the DIS data sets.
It is most common to compute the inclusive DIS cross sections using the photon's virtuality as the factorization scale, µ 2 F,DIS = Q 2 . It has been argued, however, that resummation of logarithms ln p (1/x) at x 1 improves agreement with HERA Run I+II data by several tens of units of χ 2 [31,32]. In our analysis, we observe that, by evaluating the DIS cross sections in a fixed-order calculation at NNLO accuracy, with a tuned factorization scale we achieve nearly the same quality of improvement in the description of the HERA DIS data set as in the analyses with low-x resummation [31,32]. The fit done with these modified settings is designated as CT18X. For this fit, the χ 2 (HERA I+II) reduces by more than 50 units in the kinematical region with Q > 2 GeV and x > 10 −5 , assessed in the CT18 global fit. See an illustration in Fig. 3 and its discussion in the Executive Summary II A 2.
The parametric form of the x-dependent scale µ 2 F,x is inspired by saturation arguments (see, e.g., [114,115]). The numerical coefficients in µ 2 F,x are chosen to minimize χ 2 for the HERA DIS data. At x 0.01, µ 2 F,x ≈ 0.8 Q 2 results in larger NNLO DIS cross sections than with µ 2 F = Q 2 , as it might happen due to contributions from next-to-NNLO (N3LO) and beyond. At x 0.01, µ 2 F,x numerically reduces the Q 2 -derivative of NNLO DIS cross sections. In turn, these changes result in the enhanced gluon PDF at small x and reduced gluon at x ∼ 0.01.
(ii) Unlike CT18, the CT18A analysis includes high-luminosity ATLAS 7 TeV W/Z rapidity distributions [35] that show some tension with DIS experiments and prefer a larger strangeness PDF than the DIS experiments in the small x region. Inclusion of the ATLAS 7 TeV W/Z data leads to significant deterioration in the χ 2 E values (i.e., larger S E values) for the dimuon SIDIS production data (NuTeV, CCFR), which are strongly sensitive to the strangeness PDF. One way to see this is to compare the S E distributions for the CT18 fit in Fig. 4 and the counterpart figure for CT18Z in Fig. 60 of App. A. The comparison shows that the S E values for CCFR and NuTeV dimuon data sets are elevated in the CT18Z fit, as compared to the CT18 fit, as a consequence of inclusion of the ATLAS W/Z data in the CT18Z fit. Another way to see this was carried out in Ref. [25], using the ePump program.
(iii) CT18Z represents the accumulation of these settings introduced to obtain a PDF set that is maximally different from CT18, despite achieving about the same global χ 2 /N pt as CT18. The CT18Z fit includes the 7 TeV W/Z production data of ATLAS like CT18A, but also includes the modified DIS scale choice, µ F,x , as done for CT18X. In addition to these modifications, CT18Z excludes the CDHSW extractions of the F 2 and xF 3 structure functions from νFe scattering, which otherwise would oppose the trend of CT18Z to have a softer gluon at x > 0.1, cf. Sec. II B 2. Finally, CT18Z is done by assuming a slightly higher value of the charm quark mass (1.4 GeV compared to 1.3 GeV) in order to modestly improve the fit to the vector boson production data.

28
The combination of these choices in the CT18Z analysis results in a Higgs boson production cross section via gluon fusion that is reduced by about 1% compared to the corresponding CT14 and CT18 predictions. Thus, the various choices made during the generation of four CT18(A,X,Z) fits allow us to more faithfully explore the full range of the PDF behavior at NNLO that is consistent with the available hadronic data, with implications for electroweak precision physics.

III. THEORETICAL INPUTS TO CT18
Modern global fits determine the PDFs from a large number of data points (N pt > 3600 for In the present section, we review the essential components of our theoretical setup: the goodness-of-fit function in Sec. III A, computer programs for (N)NLO computations for various processes in Sec. III B, and input parametric forms for the PDFs in Sec. III C. The explicit parametric forms for the best-fit CT18 PDFs are presented in Appendix C.

A. Goodness of fit function and the covariance matrix
The CTEQ-TEA analyses quantify the goodness-of-fit to an experimental data set E with N pt data values by means of the log-likelihood function [23] A k-th datum is typically provided as a central value D k , an uncorrelated statistical error s k,stat , and possibly an uncorrelated systematic error s k,uncor sys . Then, s k ≡ s 2 k,stat + s 2 k,uncor sys is the total uncorrelated error on the measurement D k .
T k is the corresponding theory value that depends on the PDF parameters {a 1 , a 2 , ...} ≡ a.
In addition, the k-th datum may depend on N λ correlated systematic uncertainties, and those may be fully correlated over all data points. To estimate such errors, it is common to associate each source of the correlated error with an independent random nuisance parameter λ α that is assumed to be sampled from a standard normal distribution, unless known otherwise.
The experiment does not tell us the values of λ α but it may provide the change β kα λ α of D k under a variation of λ α . Knowing β kα , one can estimate the likely values of λ α , as well as the uncertainty in the PDF parameters for a plausible range of λ α .
For those experiments E that provide β kα , we find that, at the global minimum a 0 , the best-fit χ 2 value is given as in terms of the best-fit shifted residuals, and best-fit nuisance parameters, where and These relations are derived in Appendix B.
Another instructive form expresses r i (a 0 ) in terms of the shifted data values D sh Sometimes, we take extra steps to convert the published table of correlated uncertainties into the β kα matrix formatted in accord with Eq. (1). For example, when an experiment distinguishes between positive and negative systematic variations, we average these for each data point for consistency with the normally distributed λ α . [We have verified that the choice of the averaging procedure does not significantly affect the outcomes, e.g., if a central value is shifted to be in the middle of an originally asymmetric interval, etc.] In a small number of experimental publications, only a form based on the covariance matrix (cov) ij is used in place of Eq. (1): While we can compute χ 2 directly using Eq. (8), when deriving the PDFs, we find it convenient to go back to the form consisting of the uncorrelated errors s i and the correlated contributions provided by β kα : An algorithm to construct such a representation with sufficient accuracy is presented at the end of Appendix B. In all relevant cases, we have checked that both the input covariance matrix (cov) ij and its decomposed version (9) produce close values of χ 2 . With the latter representation, we are also able to examine the shifted data values and shifted residuals, Eq. (7), to explore agreement with the individual data points.
In this article, we generally follow the CTEQ methodology and obtain r i (a 0 ) directly from the CTEQ-TEA fitting program, together with the optimal nuisance parameters λ α (a 0 ) and shifted central data values D sh i (a).

Overview
For deep-inelastic scattering observables, we perform computations using an NNLO realization [116] of the SACOT-χ heavy-quark scheme [117][118][119][120] adopted since CT10 NNLO [89]. These can be done using either the pole or MS quark masses as the input [121], grids based on the calculation presented in Ref. [122]. The impact of the NNLO contribution on the description of the charged-current dimuon DIS data is further discussed in Sec. V B 4.
The computational complexity of NNLO matrix elements precludes their direct evaluation for each fit iteration, particularly given the expansive size of the data sets fitted in CT18. Instead, for the newly included high-precision data from the LHC, ApplGrid [16] and fastNLO [15] fast tables have been generated using programs such as MCFM [42],   NNLO calculations, it was necessary to fit the K-factors with smooth curves and add a small MC error (in this analysis, equal to 0.5%) to the calculations for the high-p T Z and inclusive jet production processes, as specified below.
For the legacy data on electroweak boson production, already included in the CT14 and CT14 HERAII , we inherit the original CTEQ calculations summarized in Even with the use of stored grids for fast evaluation of the matrix elements, significant improvements on speed are needed. The CT fitting code has been upgraded to a multithreaded version with a two-layer parallelization, through a rearrangement of the minimization algorithm and via a redistribution of the data sets. As a result, the speed of calculations increased by up to a factor of 10. Details are provided in Appendix D.
We will now describe the theoretical calculations for each new LHC process included in the CT18(Z) fits.

LHC inclusive jet data
LHC inclusive jet data are available with different jet radii. We have chosen the larger of the two nominal jet radii, 0.6 for ATLAS and 0.7 for CMS, to reduce dependence on resummation/showering and hadronization effects [130]. There is a non-negligible difference at low jet transverse momentum between theory predictions at NNLO using as the momentumscale choice of either the inclusive jet p T or the leading jet p T (p T 1 ) [13] . The nominal choice adopted by the CTEQ-TEA group is to use the inclusive jet p T . We have observed that the fitted gluon PDF is not very sensitive to this choice even in the kinematic regions where the difference in NNLO predictions between these two scale choices is important. This resilience in the global fit is due both to the presence of other data in the relevant kinematic regions, and to QCD evolution.
Electroweak corrections from Ref. [131] have been applied to the jet predictions. We find that these corrections can be as large as 10% for the highest transverse momentum bin in the central rapidity region, but decrease quickly with increasing rapidity and with decreasing jet transverse momentum. Furthermore, the QCD NNLO/NLO K-factors were fitted with smooth curves, and a 0.5% theoretical error has been added to each data value to take into account the fluctuations in Monte-Carlo integration of NNLO cross sections provided by NNLOJET.

LHC electroweak gauge boson hadroproduction
The Drell-Yan theory calculations at NNLO in the CT18(Z) global analysis consist of the following: • ATLAS 7 TeV 4.6 fb −1 measurements of W ± and Z/γ * production cross sections in the e and µ decay channels [35]: the theory predictions at NLO are obtained by using APPLgrid [16] fast tables generated with MCFM [42] and validated against aMCfast [124] interfaced with MadGraph5 aMC@NLO [132]. The NNLO corrections are imported from the xFitter [133] analysis published in Ref. [35]. These corrections are obtained using the DYNNLO-1.5 code [40,41], and checked by FEWZ-3.1.b2 [45][46][47] and t MCFM-8.0 [43] codes. Some discrepancy among these codes (up to ∼ 1%) were found. However, these discrepancies do not induce significant differences in the calculated results like More details can be found in Appendix F. calculations. In addition, we have imposed the kinematic cut 45 < p ll T < 150 GeV to ensure reliability of the fixed-order calculation. The low-p ll T region is dropped due to the non-negligible contribution from QCD soft-gluon resummation, and the high p ll T region is dropped because the EW corrections are expected to be large [134,135].

Top-quark pair production
The theory predictions for top-quark pair production differential distributions at fixed order NNLO in QCD for the LHC are obtained by using fastNNLO tables [17,18]. The top-quark mass has been set to m t = 172.5 GeV. According to [125], different distributions adopt different values of the default central scale µ F = µ R = µ at which the calculation is performed. In particular, the top-quark p T spectrum is obtained using µ = 1/2 m 2 t + p 2 T,t , while the rest of the distributions are obtained with µ = 1/4 m 2 t + p 2 T,t + m 2 t + p 2 T,t . The impact of the EW corrections on tt differential distributions has been studied in [136] 35 where the difference between the additive and multiplicative approaches for combining QCD and EW corrections is investigated. EW K-factors from an analytic fit for the QCD × EW/QCD contributions are also available [137].

C. The nonperturbative parametrization form
An important source of the uncertainty in the CTEQ PDF analysis is associated with the choice of the parametric form for the fitted distributions at the lower boundary of QCD evolution, f a (x, Q = Q 0 ). The x dependence of the PDFs at this initial scale are not directly calculable in QCD perturbation theory, but are instead fitted to high-energy data after evolving from Q = Q 0 . In practice, there is limited guidance from theory as to the most appropriate parametrizations for these nonperturbative functions, and it is favorable to guarantee a maximal level of parametric flexibility without simultaneously overfitting experimental data [33]. In App. C, we present the explicit parametrization forms used in CT18.
To estimate the parametrization dependence, we repeated the CT fits multiple times using a large number (more than 100) of initial parametrization forms which have comparable numbers of fitting parameters. The final PDFs are obtained using a fixed parametrization form, but the uncertainty is computed so as to cover the bulk of the solutions obtained with the alternative parametrizations. The results of this study are illustrated in Fig. 6, which show central fits for a range of alternative fitting forms (green curves) superposed within the uncertainty band (at the 90% confidence level) for an early version of CT18.
As we increased the number of free PDF parameters, an improvement in the global χ 2 or individual S E values was typically found, so long as 30 free parameters were fitted. With more than about 30 parameters, the fits tend to destabilize, as expanded parametrizations attempt to describe statistical noise. The final PDFs are based on the parametrizations with a total of 29 free parameters. For each of the four fits, we provide twice as many Hessian error PDFs to evaluate the PDF uncertainties according to the CTEQ6 master formulas [23]. GeV. The function xf (x, Q) is plotted versus x, for flavors u, u, d, d, s = s, and g. We assume s(x, Q 0 ) =s(x, Q 0 ), since their difference is consistent with zero and has large uncertainty [55]. The plots show the central fit to the global data listed in Tables I and II, corresponding to the lowest total χ 2 for our choice of PDF parametrizations. These are displayed with error bands representing the PDF uncertainty at the 90% confidence level (C.L.).
The relative changes from CT14 HERAII NNLO to CT18 NNLO PDFs are best visualized by comparing their associated PDF uncertainties. We make a number of observations for the NNLO PDFs. We observe that the CT18 u PDF becomes slightly smaller, compared to CT14 HERAII , at almost all x values, with the largest decrease at x ∼ 10 −3 . The d PDF has increased at x ∼ 10 −3 and x ∼ 0.2, while slightly decreased at x ∼ 0.01. Theū andd distributions are both smaller at x ∼ 0.3 and larger at x ∼ 0.05, though the decrease ind is larger. Furthermore, except the d PDF at x ∼ 0.2, the error bands of u, d,ū andd are about the same as CT14 HERAII . The central strangeness (s) PDF has increased for x < 0.01 and decreased for 0.2 < x < 0.5 where the strange quark PDF is essentially unconstrained in CT18, just as in CT14 HERAII NNLO. Also, its uncertainty band is slightly larger than CT14 HERAII for x > 10 −4 , as a consequence of the more flexible parametrization and the inclusion of the LHC data. We have checked that the most important data sets that drive the above mentioned changes in the quark and antiquark PDFs are the LHCb W and Z boson data, as listed in Table II   At such high x, the CTEQ-JLab analysis (CJ15) [5] has independently determined the ratio d/u at NLO, by including the fixed-target DIS data at lower W and higher x that is excluded by the selection cut W > 3.5 GeV in CT18, and by considering higher-twist and nuclear effects that can be neglected in the kinematic range of CT18 data. Fig. 8 shows that the central prediction of CT18 differs from CJ15 at x > 0.1. Furthermore, the error band of CJ15 is much smaller than CT18 (a fact partly attributable to the ∆χ 2 = 1 criterion used in CJ15), though they are consistent. We note that here we are comparing two PDF sets at different orders in the strong coupling: CT18 is an NNLO fit, while CJ15 is an NLO fit.
Turning now to the ratios of sea quark PDFs in Fig. 9, we observe that the uncertainty Error bands of (s+ - The overall increase in the strangeness PDF at x < 0.03 and decrease ofū andd PDFs at x < 10 −3 , cf. Fig. 7, lead to a larger ratio of the strange-to-nonstrange sea quark PDFs, (s(x, Q) +s(x, Q)) / ū(x, Q) +d(x, Q) , presented in Fig. 9. At x < 10 −3 , this ratio is determined entirely by parametrization form and was found in CT10 to be consistent with x → 0, albeit with a large uncertainty. The SU(3)-symmetric asymptotic solution at x → 0 was not enforced in CT14, nor CT14 HERAII , so that this ratio at Q = 1.4 GeV is around 0.3 to 0.5 in the small x region. In CT18, we have taken a different s-PDF nonperturbative parametrization form (with one more parameter added), but still with a parametric form that ensures a stable behavior of the ratio R s , defined in Eq. (10), for x → 0, so that R s (x → 0) is about 0.7 and 1, respectively, in CT18 and CT18Z fits.
We may summarize the pulls of specific processes on the central CT18 fit as follows.
• The most noticeable overall impact of the LHC inclusive jet production on the central gluon PDF g(x, Q) is to mildly reduce it at x > 0.2 within the original PDF uncertainty band. Similar pulls are observed from jet data with decorrelation of some systematic errors, and with a 0.5% MC uncertainty on theory values. The pulls from various jet data sets on g(x, Q) neither follow a uniform trend across the whole x range nor are consistent among various measurements, as is demonstrated, e.g., by the L 2 sensitivity in Fig. 24 and LM scans in Sec. V A.
• The LHCb data, combined over all processes, have some impact on the u, d and s quarks, and pull the s(x, Q) up at small x.
• The ATLAS 8 TeV Z p T data (Exp. ID=253) have an impact on the gluon PDF at large x. Similarly to the LHC inclusive jet data, they pull the g(x, Q) down at large x.
• The ATLAS 7 TeV data on W and Z rapidity distributions (Exp. ID=248), included only in CT18A and Z, have the largest influence on the PDFs, as discussed in App. A.
The directions of the pulls are similar to LHCb.
• The LHC data on tt double differential cross sections yield a similar pull on the central fit for g(x, Q) as the LHC jet data, favoring a softer gluon at large x. However, the gluon PDF error band is mainly constrained by the inclusive jet data due to the much larger number of data points as compared to that of the tt data.

R s dependence
An observable closely connected with the flavor-structure of the nucleon sea is the strangeness suppression factor, R s measures the x and Q dependence of the breaking of flavor-SU(3) symmetry, with older analyses typically fixing R s = 0.5. More recently, a number of previous CTEQ studies [55,138] examined then-current constraints on the behavior of R s , particularly as driven by the neutrino-induced SIDIS dimuon production measurements, done by the CCFR and NuTeV Collaborations. These works found significant evidence of an independent x dependence for , distinct fromū+d, but were unable to exclude a vanishing strangeness momentum fraction asymmetry, In the present work, we do not separately parametrize s(x) ands(x), but investigate the impact of the newlyfitted LHC data on the total strange distribution, s + (x), and flavor asymmetry ratio, R s (x) in the presence of the older neutrino dimuon data.
Among these experiments, we observe an especially marked sensitivity to R s from the LHCb and ATLAS W /Z production data (for the latter data set, examined with the CT18A and CT18Z fits). The cumulative effect of these experiments is to enhance the magnitude of R s at small x 0.03 for Q ∼ 1.4 GeV, as can be seen in the lower panels of Fig. 9; this effect extends to a corresponding increase in R s at the common reference of x = 0.023, Q 2 = 1.9 GeV 2 . We find this result to be stable against the specific parametrization choice, and is generally characterized by the relative increase of s =s and the comparative constancy of u,d at low x in the update to CT18. At high x, on the other hand, we find a more pronounced influence of the parametrization choice on the resulting behavior of R s , as indicated by the numerical instability we obtain in LM scans for R s at high x, as shown in     The scan exercise can also be carried out at NLO in α s , as we show in Fig. 13. In fact, any difference between the NLO and NNLO results can serve as a partial estimate of the theoretical uncertainty of its determination. Although the uncertainty is similar to that obtained at NNLO, the central value is slightly higher: α s (M Z , NLO) = 0.1187 ± 0.0027.
We note that the qualitative interplay of the experiments with leading sensitivity to α s (M Z ) is much the same at NLO as found at NNLO, with the combined HERA (Exp. ID=160) and BCDMS F p 2 data (Exp. ID=101) again preferring lower values, while the ATLAS 7 TeV jet data (Exp. ID=544) and 8 TeV Z p T data (Exp. ID=253) pulling in the opposing direction, but more strongly at NLO than at NNLO.
To summarize, we find that the CT18 data set prefers a larger value of α s (M Z ) and a marginally smaller nominal uncertainty than in CT14. is presented in Fig. 75 in App. A. higher.
An illustration of the observed trends can be viewed in Fig. 14 for the gluon-gluon luminosity, this suppression can be as large as ∼ 4% between 100 GeV and 1 TeV.
In Figs. 16 and 17, we illustrate a comparison between parton luminosities obtained using CT18 and PDFs from other groups such as CJ15 [5], MMHT14 [2], and NNPDF3.1 [3], at the LHC collider energy of 14 TeV. For such a comparison, we adopt the re-scaled 68% C.L. for the CT18 PDFs to match the convention of the other groups. In particular, the comparison shown in Fig. 16  Knowledge of the integrated PDF Mellin moments has long been of interest, both for their phenomenological utility, and for their relevance to lattice QCD computations of hadronic structure [21,140]. In the case of the former, PDF moments can serve as valuable benchmarks for the purpose of comparing various global analyses and theoretical approaches, and can also be informative descriptors of the PDFs themselves. This follows especially from the fact that numerical results obtained for PDFs of a given order are connected with the x dependence of the underlying parton distribution, with, in general, higher-order Mellin moments mostly determined by the PDFs' high-x behavior.
Integrated moments can in general be evaluated for practically any phenomenological PDF from its underlying distribution, provided the moment in question is convergent over the full range of support. However, in this analysis we concentrate special attention on with n = 1 for the gluon, as well as for the quark distributions, where we denote the charge conjugation-even (odd) quark combinations as q ± = q ± q. We primarily consider these specific PDF moments of Eq. (11) with n = 1 and Eq. (12) for compatibility with lattice QCD determinations, which are only able to compute the odd (n = 1, 3, . . . ) moments of q+q-type distributions and even (n = 2, 4, . . . ) moments for q −q-type distributions. This follows from the fact that lattice calculations extract the integrated Mellin moments from hadronic matrix elements as, In Eq.   Fig. 19. The CJ15 uncertainty for Observations. In general, we observe concordance among the moments of the light distributions, including those of the isovector (i.e., u − d) combination, x 1,3 u + −d + and Notably, the CT results for the first isovector moment, x u + −d + ∼ 0.156−0.159, are marginally larger than those obtained under the other fits considered here, which produce 152, but are nevertheless in close agreement at the 1σ level. Similarly, we recover very robust agreement for the first moment of the gluon PDF, which can be understood to carry ∼ 41% of the proton's longitudinal momentum at the scale Q = 2 GeV.
We find a slightly smaller total contribution to the momentum sum rule from the gluon under CT18Z NNLO, which results in x g = 0.407 (8), in agreement with the default CT18 NNLO calculation, x g = 0.414 (8). This is consistent with the modest reduction in the central gluon shown for CT18Z in the lower-right panel of Fig. 7.   [141] and 1994 [142], respectively. These were based on direct quark-parton model extractions of the flavor asymmetry PDF from the deuteron-to-proton structure function ratio measured at Q 2 = 4 GeV 2 for a range of x 0.7. While all but the highest x bin in this data set is consistent with CTEQ kinematical cuts, the very low Q is exactly at the boundary of the Q cut, and likely subject to substantial higher twist and mass corrections, especially for the higher x bins.
longitudinal momentum, we again find in general strong convergence among our new global analysis and the results of previous and other fits. This is especially true for the total u + and d + first moments, for which we find concordance at x u + ∼ 0.35 and x d + ∼ 0.193−194. The situation is similar for the total nucleon strangeness momentum, but with a somewhat greater quantitative spread about x s + = 3 − 4%. For CT18 NNLO, we obtain x s + = 3.3 ± 0.9% -similar to CT14 HERAII . In shifting to CT18Z, a 24% larger nucleon strange content is preferred, but with comparable error.
In addition to the first moments of the quark and gluon densities, x q + ,g , we also evaluate the second moments of select q −q quark asymmetries according to Eq. (12), finding for x 2 u − ,d − very close alignment among CT18(Z) and previous calculations. Recent CT fits and CJ15 do not independently parametrize s vs. s, and we therefore omit entries in Table VI for x 2 s − . Results on the integrated PDF moments are also of interest to phenomenological sum rules -for instance, the Gottfried Sum Rule [143], which relates x −1 -weighted moment of to flavor-symmetry violation in the light quark sea via the breaking of the SU (2)  analyses [141,142]; although the uncertainties of the CT determinations are considerably larger, as indicated by the black, red, and green bars in Fig. 19; the two nested bands correspond to the extracted values of 1 d −ū found in Refs. [141,142]. We note that CT uses a significantly more flexible parametrization for the light-quark sea (with 11 parameters for the combinedū andd-PDFs in CT, compared with 5 parameters in MMHT14 ford −ū [2] and 5 parameters in CJ15 ford/ū [5]), with no constraint on the sign ofd−ū, as can be deduced from thed/ū(x, Q = 1.4 GeV) ratio plot shown in the upper-left panel of Fig. 9. We therefore find that, with this flexibility, especially in the low-x region important for 1 d −ū , modern high-energy data still allow a broad range for the zeroth moment.
We may extend the analysis of the d = u breaking to the SU(3) sector, by analyzing the ratio of the first moments of the distributions appearing in Eq. (10) leading to the strange suppression factor moment ratio, as illustrated in Fig. 18. We detail the numerical results for this quantity in the final row of  Refs. [144,147,148], while the corresponding result for x 2 u − −d − is an average over the result in Ref. [149] and two separate calculations reported in Ref. [150]. strangeness suppression scenario, κ s = 0.5. In moving from CT14 HERAII to CT18, there is a modest enhancement in the preferred central value and related growth of the associated uncertainty, which shifts from κ s CT14 HERAII = 0.46±0.13 to κ s CT18 = 0.49±0.16, in very close agreement with MMHT14, in particular. The removal of CDHSW and inclusion of the AT-LAS W, Z production data, as well as other changes leading to CT18Z result in a significant increase in the ratio, which increases to κ s CT18Z = 0.61±0.14, with a small contraction of the PDF uncertainty, as compared to CT18. The growth in the value we find for κ s for CT18Z is evidently driven by an enhancement in the strangeness sector, given the corresponding enlargement of x s + described above. This increase in the uncertainty brings it in the same direction as NNPDF3.1, which also prefers a slightly larger value for κ s .
Recently, a first lattice calculation of κ s was reported by the χQCD collaboration in Ref. [146], which found κ s (Q = 2 GeV) = 0.795 ± 0.079 (stat) ± 0.053 (sys). Indeed, while this result lies just beyond the upper periphery of the values preferred by typical phenomenological fits, κ s ∼ 0.5, it agrees at the 1σ level with the CT18Z result that follows the inclusion of 7 TeV inclusive W, Z production data taken by ATLAS.
Examining the entries appearing in the rightmost column of Table VII, it is evident that PDF moments determined on the QCD lattice have a general tendency to overestimate the values extracted phenomenologically. In some circumstances -e.g., for the isovector u−d moments, or for the total d-quark momentum, x d + -the lattice uncertainty is sufficiently large that there is nevertheless agreement with global analyses. In contrast, the moment associated with the gluon momentum x g can also be independently computed on the lattice, with results for this quantity often substantially smaller than those determined from PDF fits; this includes the calculation of Ref. [144], which still recovered the momentum sum rule, x u + +d + +s + + x g = 1, up to uncertainties by merit of the larger total quark momenta.
Schematically, the PDF moments are extracted on the lattice from the ratio of 3-point to 2-point correlation functions [140,151]: where the B a,b are baryon interpolating operators, t the source-sink Euclidean time separation, and τ the Euclidean time associated with the operatorÔ insertion noted in Eq. (13).
For the lower moments of the nucleon parton distributions, the lattice output is substantially governed by the interplay between excited-state contamination of the correlation functions, which in general depend on Euclidean time as ∼ exp (−m i t), and the lattice signal-to-noise ratio, which goes as S/N ∼ exp (−(E N − [3/2]m π )t). As such, lattice calculations at physical pion mass (or chiral extrapolations thereto) lead to more rapid deterioration of the signal-tonoise at precisely the larger lattice times at which contributions from nucleon excited states are relatively suppressed. The subtle relationship between these lattice effects (in addition to other systematic artifacts) complicate any straightforward interpretation of the presently large(small) lattice results for the quark(gluon) PDF moments shown in Table VII. 60

V. QUALITY OF THE FIT TO DATA
The CT18 global analysis includes a wide range of data from Run-1 of the LHC, in addition to the extensive collection of data used in the previous CT14 analysis with the combined HERA measurements. This expansive data set was introduced in Sec. II B and   Tables I-II,  Then, we minimize χ 2 with respect to the parameters a of the functional forms of the parton distribution functions. We arrive at the best-fit χ 2 given by Eq. (2) as the sum of (D sh i (a 0 ) − T i (a 0 )) 2 /s 2 i and squares of optimal individual nuisance parameters λ(a 0 ). Here T i is the i-th theory prediction, D sh i denotes the respective data value shifted by the optimal systematic displacements of the nuisance parameters; s i is the published estimate for the total uncorrelated error.  Tables  I and II. [34].
If all deviations of theory from data are purely due to random fluctuations, one would expect to recover an empirical distribution of S E that is close to the standard normal distribution N (0, 1) (with mean 0 and standard deviation 1) for any N pt,E . In fact, any recent global fit renders an S E distribution that is statistically incompatible with N (0, 1) [33], indicating that many experiments are underfitted or overfitted.
For the CT18 NNLO fit, the observed S E distribution shown in Fig. 4 is most compatible with N (0.6, 1.9). The probability that is compatible with N (0, 1) is very small (p = 2.5·10 −5 according to the Anderson-Darling test [33]). In the figure, we labeled the experiments with the largest deviations from S E = 0. These are the combined HERAI+II data set on inclusive DIS [27] with S E ≈ 5.7, which provides the dominant constraints on the PDFs and must be retained in the global analysis despite the quality-of-fit issues discussed in Sec. II A 2; and the CCFR measurement [70] of the structure function xF 3 (x, Q) in charged-current DIS on iron, which has an unusually low χ 2 /N pt ≈ 0.4 for the central fit, but does constrain the PDF uncertainty for some flavors, as can be seen, e.g., in the LM scans presented in the next section.
We also note that the newly LHC Run-1 data sets, indicated by the light green color in  Table II. Two squares and two stars indicate the S E values for the NuTeV dimuon and CCFR dimuon data, respectively, which we highlight for special attention given the importance of these data for probing the strangeness PDF. An analogous plot for the alternative CT18Z fit in Fig. 60 shows increased S E values for the CCFR and NuTeV experiments, as compared to the CT18 fit, because of the conflicting pull of the ATLAS 7 TeV W/Z production data.   At the global minimum itself, the χ 2 E for such an experiment may be elevated by up to tens of units.
In Fig. 20, for instance, we show two LM scans associated with the gluon density, g(x, Q).
In the left panel, the LM scan probes the pulls of the most sensitive measurements to the Higgs-region gluon PDF, which contributes to Higgs boson production through the predominant gg → H channel, especially in the neighborhood of x = m H /(14 TeV) ∼ 0.01 and for Q ∼ m H . Evidently, most constraints arise due to HERA inclusive DIS data as well as the LHC jet data. On the right plot for x = 0.3, the constraints spread over more data sets, e.g., from top-quark production data, Z boson p T data, and some fixed target DIS data.
For strange quark density, which is also shown in Fig. 22, the only significant constraints arise from the NuTeV and CCFR dimuon data, and also the ATLAS 7 TeV W and Z data if included. That is true for both small and large-x values also because of the limited parametrization form of strangeness. For other sea quarks and alsod/ū ratios, the LHCb data and CMS W boson charge asymmetry data play important role at small-x, as can be seen from Fig. 23. On another hand, at large-x the flavor separation mostly rely on E605, E866 and NMC deuteron data.
The advantage of this approach is its systematic and robust nature. On the other hand, this calculation requires dedicated refits of the PDF parameters subject to the LM constraints for a sufficiently large collection of LMs to map the nontrivial dependence of the experiments' χ 2 on the local PDF behavior -a fact which makes the LM scans computationally expensive.

L 2 sensitivity analysis
A technique complementary to the LM scans explored in Sec. V A 1 is the calculation of the L 2 sensitivity. The L 2 sensitivity was first introduced in Ref. [21] for the purpose of analyzing the interplay among the pulls of the CT18(Z) data upon the fitted PDFs.
While the LM scans offer the most robust approach for exploring possible tensions among fitted data sets in a given analysis, they are very computationally costly to evaluate, and generated at specific choices of x and Q. As we explain here, the L 2 sensitivity can be rapidly computed and provides a strong approximation to the ∆χ 2 in a particular global analysis. Moreover, the L 2 sensitivity can be readily computed across a wide range of x, allowing the ∆χ 2 variations shown in the LM scans to be visualized and interpreted in an x-dependent fashion. We stress that the qualitative conclusions revealed by consideration of the L 2 sensitivities discussed and presented below are consistent with the picture based on the LM scans themselves. Although the L 2 sensitivities may not always provide the same numerical ordering as the LM scans for the sub-dominant experiments, it offers a complementary information over broader reaches of x not entirely captured by individual LM scans.
We work in the Hessian formalism [22,23,59] and compute the L 2 sensitivity S f,L2 (E) for each experiment, E, as which yields the variation of the log-likelihood function χ 2 E due to a unit-length displacement of the fitted PDF parameters away from the global minimum a 0 of χ 2 ( a) in the direction of ∇f . The PDF parameters a are normalized so that a unit displacement from the best fit in any direction corresponds to the default confidence level of the Hessian error set (90% for CT18).
This displacement increases the PDF f (x, Q) by its Hessian PDF error ∆f , and, to the extent its PDF variation is correlated with that of f (x, Q) through the correlation angle In order to understand the interplay and competing pulls among the fitted CT18 data, we compute the L 2 sensitivity for several representative examples of importance to LHC phenomenology. A full collection of the L 2 sensitivity plots for CT18(Z) PDFs and PDF ratios can be viewed at [152]. Analogous calculations are shown for the alternative CT18Z fit in Sec. A 4.
In Fig. 24, we show the L 2 sensitivity of the CT18 data on the gluon PDF at fixed Q = 100 GeV, plotting curves for those experiments that satisfy |S f,L2 (E)| ≥ 4 for any value  Fig. 20 (left). Of the newly fitted LHC Run-1 data in CT18, the 7 TeV inclusive jet production data taken by ATLAS (Exp. ID=544) shows strongest overall pull at x = 0.01, S g,L2 (E) ∼ 5, although the corresponding pull of the CMS 8 TeV jet data (Exp. ID=545) is still larger for slightly higher x ∼ 0.05 − 0.1. The effect of the Run-1 jet data is nonetheless overshadowed at x = 0.01 by the pull of the combined HERA Run 1+2 data (Exp. ID=160), for which we find S g,L2 (Exp. ID = 160) ≈ 7 for Q = 100 GeV.
These observations are consistent with our findings based on the LM scans appearing in Sec. V A 1, as typified by Fig. 20 (left panel), wherein we identified the same experiments as imposing the most stringent constraints upon g(x = 0.01, Q = m H ). Notably, Fig. 24 also makes evident the competing tensions among the fitted data sets in the Higgs region.
While, for the most part, the data sets fitted in CT18 pull g(x = 0.01, Q = m H ) in the same direction as the LHC jet production and HERA data mentioned above, a few experiments exhibit competing pulls, with PDF variations in their χ 2 E being anti-correlated with the fluctuations in the Higgs-region gluon PDF. The most pronounced of these are the ATLAS 8 TeV Z p T data (Exp. ID=253) and the fixed-target pp Drell-Yan cross section information recorded by E866 (Exp. ID=204), for both of which we find S g,L2 ≈ −4 under CT18.
The pulls and tensions on the behavior of the high-x gluon distribution can be similarly understood by examining Fig. 24 at x 0.3, which is also marked by a number of very strong competing pulls. Among the more striking of these is the opposing alignment of the CMS 7 and 8 TeV jet-production data.

B. Description of data sets fitted in CT18
In this subsection, we illustrate the ability of CT18 to describe the individual experiments included in this analysis, with particular attention paid to the newly included LHC Run-1 data. We organize this discussion according to the specific physical process.

Vector boson production data
Legacy Tevatron data. The charge asymmetry data at D0 Tevatron Run-II are described well within the CT18 fit. Fig. 25 shows a data versus theory comparison for the electron charge asymmetry as a function of the absolute value of the electron pseudorapidity. TeV is in Fig. 28  improved after dropping the first rapidity bin data.
The general quality of the CT18 fit to these data can be quantified by the residual plots, as shown in Fig. 29. For example, the third panel illustrates the distribution of values of LHCb 8 TeV W ± and Z data (Exp. ID=250). It indicates that there are a few data points with large values in the Exp. ID=250 data set. As expected, the large residuals result from the first few rapidity bins near y = 2 in the W ± and Z data, and from 3.5 y Z 4 between 3.5 and 4 in the Z data. Another useful criterion is the examination of the distribution of nuisance parameters needed to fit the Exp. ID=250 data, which is shown in Fig. 30. As show in Fig. 31, the distribution of nuisance parameters deviates from the normal distribution. In particular, there are two nuisance parameters with large values (±3) which represent the pulls of these data on the CT18 PDFs including u, d, s,d, andū, at Q = 100 GeV. The theoretical predictions, compared with the shifted and unshifted data are shown in Fig. 32. We see that all the experimental data are fitted well within the 68% C.L. PDF uncertainty. The corresponding residuals and nuisance parameters are shown in Fig. 33.
Both are reasonably compatible with the normal distribution with mean 0 and standard deviation 1.
In the CT18(Z) analysis, we have also included the transverse momentum (p T ) distributions of lepton pairs produced in Z decays measured at ATLAS at a center of mass of 8 TeV   TeV W -lepton charge asymmetry data (Exp. ID=249).
(Exp. ID=253). The theoretical predictions for these data are obtained based on the NNLO fixed-order calculations for Z+jet production. We want to emphasize that we have imposed a kinematic cut 45 < p Z T < 150 GeV to remove both low-and high-p T data. The low-p T data are dropped because of the missing resummation effects in our fixed-order calculation.
The high-p T data are cut off because the EW corrections are non-negligible, which will be discussed in Sec. V C.
We generated the NLO APPLgrid files with MCFM. The perturbative K-factors for these data are computed directly as the ratio of the NNLO and NLO cross sections as obtained from Refs. [50,51,[153][154][155][156]. To account for non-negligible fluctuations in the NNLO theoretical prediction, we have included an additional 0.5% theoretical Monte-Carlo uncertainty, The CT18 NNLO theoretical predictions are compared to the ATLAS 8 TeV data in Fig. 34.
We see that for CT18 NNLO to describe the data in all 3 invariant-mass bins, substantial systematic shifts are required. The distributions of the residuals and nuisance parameters of this data are shown in Fig. 35. There is an overall systematic shift when comparing the unshifted data to the theoretical predictions in Fig. 34. The corresponding nuisance parameter associated with this shift, due to the luminosity error, is λ lumi = 2.02, as can be seen in the upper tail of the distribution plotted in the right panel of Fig. 35. We have investigated the scale dependence by varying the renormalization and factorization scales independently by a factor of 2, while restricting their ratio to be We find that the optimal scale choice to describe the unshifted data is given by (see Sec. V C), pulling the theory further away from the data. Instead, the systematic shift in the normalization can possibly be ascribed to the missing higher-order (N 3 LO) corrections, implied by two observations. First, the NNLO corrections to the Z p T are generally as large as 10%, which indicates slow convergence of the perturbative expansion. Second, the large scale uncertainty (about 3-4%) is also an indication that the missing higher-order effects may be significant.

Jet data
In addition to the data from the Tevatron Run-II fitted in CT14, CT18 now incorporates inclusive jet production data from Run-1 of the LHC, as described in Sec. II B 4. Historically, inclusive jet production has played an important role in constraining the gluon density, g(x, Q), in previous phenomenological analyses. This fact was evidenced by the impact that the older Tevatron jet data had in the previous CT14 global analysis. CT18 now implements higher integrated luminosity inclusive jet production data measured by the ATLAS and CMS collaborations at the LHC.
Tevatron Run-II data. First, we examine the fits to the Tevatron Run-II jet data.
The comparison shown in Fig. 37 indicates that the CDF Run-II high p T jet data prefer a harder gluon at large-x provided by the CT18 PDFs, while the D0 Run-II jet data, depicted in Fig. 38, show better agreement with the new LHC jet data included in the CT18 fit, and favor a softer gluon at large-x values as compared to the CDF Run-II jet data.
Run-1 LHC data. The CT18 fit can describe the LHC CMS and ATLAS jet data, after including a careful treatment of the correlated systematic errors, using the decorrelation prescriptions laid out in App. E for the ATLAS data, as well as the inclusion of a 0.5% overall uncorrelated systematic error for all the LHC jet data, as discussed in Sec. II A number of interesting patterns emerge for each data set. For example, in the ATLAS inclusive jet data at 7 TeV, the best fit requires the correlated errors to shift the raw data downward in the smaller rapidity regions, but to shift upward the raw data for high rapidities.
For the 7 inclusive jet production data from ATLAS, as well as the corresponding 7 TeV jet data from CMS, we track the values of the nuisance parameters λ α for k ∈ [1, N λ ], where N λ represents the number of correlated systematic uncertainties. The fitted values of the nuisance parameters obtained with CT18 for each of the three inclusive jet production data sets are plotted in the histograms of Fig. 40. In general, the values of these nuisance parameters are distributed narrowly about |λ α | ∼ 0. However, the values of the nuisance parameters are substantially larger for some correlated systematic error sources, which we describe in greater detail in App. E.

Top-quark pair production data
The two tt data sets included in the CT18 global analysis are well-described, as shown by the values of χ 2 and effective Gaussian variable S E given in the latter two rows of   Tables 5 and 7, respectively, of Ref. [98]. We have checked that including the covariance matrix of statistical correlated uncertainties (listed in Table 6  In Fig. 44, the top-quark p T distribution at CMS is fitted reasonably well across the four rapidity bins examined here. We find modest deviations between theoretical predictions and the (un)shifted data for some points in the intermediate rapidity bins, 0.75 < |y t | < 0.85 and 0.85 < |y t | < 1.45, contributing to the somewhat broader distribution of residuals for this experiment plotted in the upper-left panel of Fig. 46. Notably, the effect of correlations in fitting the CMS data is relatively minimal, given the fact that shifting these data based on nuisance parameters is generally small as evidenced by the similarity between the black (unshifted) and blue (shifted) data in Fig. 44 In contrast, achieving a very good description of the analogous ATLAS p T,t and m tt distributions shown in Fig. 45 critically depends on the use of nuisance parameters to compensate for correlated systematics as seen in Fig. 45.

Dimuon production
For CT18, we have closely examined theoretical predictions on cross sections of various processes that are sensitive to the strange-quark PDFs, e.g., charm-quark production cross sections in neutrino deep-inelastic scattering. In the CT14 NNLO analyses, we used the charm-quark production cross section calculated at NLO in QCD [158][159][160]   NNLO calculation into our public release of CT18 NNLO PDFs for the following reasons.
A consistent treatment of charm-quark mass effects at NNLO requires non-trivial matching of the charged-current results from FFN calculations to an ACOT-like VFN scheme, which is not available yet. Furthermore the charm-quark production cross sections as measured by CCFR and NuTeV are extracted from dimuon cross sections using a MC generator with NLO precision only. However, the discussion in Ref. [161] shows that, for the kinematics involved in CCFR and NuTeV dimuon experiments, the differences at NNLO between the results from the FFN scheme and of any VFN schemes are expected to be much smaller as compared to the precision of the experimental data. Thus, in the CT18 analyses, we perform alternative fits by simply using the FFN calculations of charm-quark production cross sections at NNLO from Ref. [161].
We compare CT18 and CT18Z NNLO fits, both using the charm-quark production cross sections at NLO, to their alternatives using the NNLO cross sections mentioned above. In the case of CT18, the global χ 2 is reduced by 6 units, while the reduction in the alternative fit is of a couple of units for the total χ 2 of dimuon data. . While for CT18Z, the global χ 2 and the total χ 2 of dimuon data in the alternative fit are reduced by 11 and 8 units, respectively. In both cases the NNLO predictions do provide slightly better agreement with the data though the improvement is mild, especially considering the large tolerance used in CT18(Z) analyses.
We present the strange-quark PDF and the strange to non-strange sea-quark ratio in Fig. 47, comparing the nominal fits and their alternatives using NNLO predictions. We observe an increase of the strange-quark PDF at x close to 0.1 in the CT18 alternative fit, similar to the PDF profiling results in Ref. [161], as understood due to the negative NNLO QCD corrections in the same x region. On the other hand, the changes are small compared to the PDF uncertainties. This can also be inferred from the relative stability of the χ 2 values. In the CT18Z fit, the impact of NNLO corrections on the PDFs is even smaller as compared to the CT18 fit, because the strange-quark distributions are mostly driven by the ATLAS 7 TeV W/Z data.
The NNLO corrections do push the strangeness towards the direction preferred by the ATLAS 7 TeV W/Z data. This can be clearly seen by extending the ePump study done in Ref. [25], in which it was shown that the s(x, Q) PDF in the CT14 HERAII fit is mostly constrained by the NuTeV and CCFR dimuon data, cf. Fig. 31 of [25]. In Fig. 48, we compare the strange-quark PDF obtained from four different fits. The base PDF set (labeled as "CT14HERA2mDimu.54") is the global fit with NuTeV and CCFR dimuon data removed from the CT14 HERAII data set. Adding back those four dimuon data sets, with NLO and NNLO predictions, yields the "CT14mDimupDimuNLO" and "CT14mDimupDimuNNLO" fits, respectively. The "CT14mDimupDimuATL7ZW" fit is obtained by adding the ATL7ZW data set to the "CT14HERA2mDimu.54" fit. While the NNLO prediction yields an s PDF closer to that constrained by ATL7ZW, the improvement is still too weak to resolve the tension between the ATL7ZW and dimuon data sets. This conclusion holds in the CT18 and CT18Z fits.

C. Electroweak corrections
In this subsection, we present a description summary of the electroweak (EW) corrections to the LHC data with upper bounds on them listed in Table VIII. The largest EW correction (8%) to the inclusive jet data occurs in the highest p T jet bins in the central jet rapidity region. We have already taken it into account in the CT18 fitting, using the multiplicative K-factor method. The EW corrections to the tt production are already investigated in Sec.   high-mass Drell-Yan M ll ∼ 1 TeV +5%(PI)-3%(EW) FEWZ Z p T p T ∼ 1 TeV -30% [135] III. B, which shows that the largest EW corrections occur to the p T (t) distribution. At the largest p T (t) data about 500 GeV, the EW correction is around -5%, it dies out rapidly with the decreasing p T (t). For other observables, the EW corrections can be negligible, as compared to the experimental uncertainty.

89
The EW corrections to the inclusive W + , W − and Z/γ * production data have been investigated in Ref. [35] using the MCSANC framework [162]. For W + and W − production, the EW corrections were found to be −0.4% and −0.3%, respectively. For the neutral-current for neutral current DY cross sections. Due to this reason, we decided not to include the lowand high-mass and forward Z-peak DY data in the CT18(Z) fits.
We did not include the ATLAS 8 TeV high-mass Drell-Yan data [163] in our CT18(Z) fitting, due to the non-negligible EW corrections and photon-induced (PI) contribution.
Their effect can be calculated using the FEWZ program, and is shown in Fig. 49. We find that, at a large invariant mass (M ll ∼ 1 TeV), the PI contribution can be as large as 5%, computed with LUXqed17 plus PDF4LHC15 [8]. In comparison, the EW corrections are about 3%, with a negative sign. The partial cancellation of the PI contribution and EW corrections yields an increase in the cross section by less than 2%. With the ePump program, we have also checked that the impact of these data on the CT18 fits is very mild and does not alter the resulting PDFs.
The only Z p T distribution included in CT18(Z) fits is coming from the ATLAS 8 TeV measurements. We have dropped the high-p T data by imposing a kinematic cut p Z T < 150 GeV, because the missing EW corrections to the high-p T data are non-negligible. In general, the EW corrections are negative. In terms of Refs. [134,135], the NLO EW corrections can be as large as a few tens of percent, when p Z T M Z , due to the electroweak Sudakov logarithms. In the fitted region of p Z T , between 45 GeV and 150 GeV, the EW corrections are found to reduce the cross sections by several percent, thereby pulling the theory predictions further away from the ATLAS 8 TeV data.  The degree of anti-correlation found in the ggH and Z boson cross sections decreases as the LHC energy increases. On the other hand, Fig. 51 shows that the electroweak gauge boson cross sections are highly correlated with each other at the LHC. Generally speaking, the prediction of CT18 is closer to CT14 HERAII , and the largest difference occurs between CT18Z and CT14. Furthermore, the CT18X prediction is closer to CT18Z for the electroweak gauge boson productions, cf. Fig. 51, but not for the ggH or tt inclusive cross sections, cf. Fig. 50.    As described in the previous sections, we have performed a global fit of the experimental data to the NNLO predictions, which results in the CT18 and its family PDF sets. This also applies to the vector boson differential cross sections measured at the Tevatron and the LHC. In contrast, we had compared this type of precision data to the ResBos predictions, which include effects from multi-gluon emission [49], to produce the CTEQ6.6, CT10 and CT14 PDFs. Hence, it is interesting to compare the vector boson differential cross section measurements to the predictions of CT18 (and CT18Z) PDFs with NNLO and ResBos calculations. As an example, we compare to the LHCb 7 TeV W and Z boson differential distributions [91]. As shown in Fig. 52, both NNLO and ResBos predictions agree within the experimental error bands.
C. W plus charm-jet production at the LHC TeV W + +c-jet (left) and W − +c-jet (right) production, respectively.
As concluded in the previous section, the s-PDF of CT18 differs from CT14 in the small- Indeed, we observe a constant upward shift in the predictions based upon CT18Z compared with CT18 in both panels of Fig. 53 for W + and W − . Interestingly, for W + production, the CT14 HERAII predictions tend to overshoot those of CT18(Z) in the large-rapidity region, a feature which agrees with the large-x behavior of the strangeness distribution seen in the lower-left panel of Fig. 7. On the other hand, for W − production in the right panel of Fig. 53, the CT14 HERAII predictions lie well below the CT18(Z) predictions over the full plotted range due to an enhanced d-quark distribution (which contributes via a Cabibbo-suppressed channel) of CT18(Z) as preferred by the LHCb 7 TeV W/Z data as shown in the upper-right panel of Fig. 22. Fig. 54 also shows that these data are most correlated with the s-PDF in the x-range from 0.01 to 0.2, where CT18 and CT18Z differ from CT14 HERAII , cf. Fig. 7.
D. Top quark pair differential distributions at the 13 TeV LHC In this section we present the data-to-theory comparison for the CMS 13 TeV measurements of top-quark pair production differential cross sections in the dilepton channel [105].
The QCD theoretical predictions at NNLO in QCD are obtained by using fastNNLO ta- the slope between the theory and unshifted experimental data for both dσ/dp t T and dσ/dm tt . Those differences are accommodated by systematic shifts of the data, resulting in a good χ 2 after all uncertainties are taken into account. We notice that, in the case of the p T spectrum, the theory prediction obtained with CT18Z NNLO gives a slightly better description of the data at large p T .
The impact of the electroweak corrections is illustrated in Fig. 58. These corrections are included as K-factors using the multiplicative scheme according to Ref. [136]. They were downloaded from Ref. [137]. Large EW effects show up in the high p t T tails. However in the p T range 1 -500 GeV shown in the figures, the EW corrections are not larger than 3-4% in most cases. If one considers higher-p T regions, K-factors would be much larger there.
The EW corrections minimally improve the agreement of theory and data for the top-quark p T and m tt distributions. The χ 2 /d .o.f . of the NNLO QCD + NLO EW prediction using CT18 PDFs agrees well with the values presented in Table 49 of Ref. [105]. For all other distributions, the EW corrections are negligible for the kinematic range studied.
Among various one-dimensional tt differential distributions, the distribution of the topquark pair rapidity, y tt , shows a good agreement between the CMS data and CT18 predic- tions. To examine how this data could modify the CT18(Z) gluon PDFs, we use the ePump program [25] to update the CT18(Z) PDFs, after including the CMS 13 TeV y tt data in the fit. As shown in Fig. 59, the updated gluon-PDF error band (labeled as CT18pCMS13ytt) is very slightly reduced for x from 0.1 to 0.4 in both cases. Further discussion about these data sets will be presented elsewhere.

VII. DISCUSSION AND CONCLUSIONS
In this paper, we have presented the CT18 family of parton distribution functions (PDFs), including the CT18Z, CT18A, and CT18X alternative fits. CT18 is the next generation of NNLO (as well as NLO) PDFs of the proton from a global analysis by the CTEQ-TEA group.
This represents the next update following the release of the CT14 and CT14 HERAII NNLO distributions, the latter of which was prompted by the release of precision HERA I and II combined data after the publication of CT14. Although some of the early 7 and 8 TeV LHC Run-1 data, including measurements of inclusive production of vector bosons [62,[83][84][85] and jets [63,64], were included as input for the CT14 fits, CT18 represents the first CT analysis that substantially includes the most important experimental data from the full Run-1 of the LHC, including measurements of inclusive production of vector bosons, jets, and top quark pairs at 7 and 8 TeV. Detailed information about the specific data sets included in the CT18 global analysis can be found in Tables I and II,  The challenge from the experimental side is to select and implement relevant and consistent data sets in the global analysis. Specifically, we have included processes that have a sensitivity for the PDFs of interest, and for which NNLO predictions are available. For example, we include as large a rapidity interval for the ATLAS jet data as we can, using the ATLAS decorrelation model, rather than using a single rapidity interval. We noted that using a single rapidity interval may result in selection bias. The result may be a larger value of χ 2 /N pt , due to remaining tensions in the ATLAS jet data, as well as reduced PDF sensitivity compared to the CMS jet data, cf. Sec. II.B.2. Similarly, to incorporate the tt differential cross section measurements into the CT18 global analysis, we use two tt single differential observables from ATLAS (using statistical correlations) and doubly-differential measurements from CMS in order to include as much information as possible. Again, there is a risk of bias if using only a single differential distribution, as some of the observables  Tables IV and V. For example, a nonnegligible difference was found at low-jet transverse momentum between theory predictions at NNLO using as the momentum-scale choice either the inclusive jet or the leading-jet transverse momentum [13]. The nominal choice adopted by the CTEQ-TEA group is to use the inclusive-jet p T . We have observed that the fitted gluon PDF is not very sensitive to The PDFSense program [20] provides an easy way to visualize the potential impact of data on PDFs in the x and Q plane. In addition, a simple L 2 sensitivity variable [21]  The inclusive of new data and theoretical advances have resulted in the following changes in CT18, as compared to CT14: 1) a smaller g(x, Q) for x ∼ 0.3 (mainly due to ATLAS and CMS 7 TeV jet data and ATLAS 8 TeV Z p T data, with some tension found between CMS 7 and 8 TeV jet data); 2) some changes in u, d,ū andd at small x, such as a larger d and d/u and a smallerd/ū for x ∼ 0.2 (mainly due to LHCb W and Z rapidity data and CMS 8 TeV W lepton charge asymmetry data); and 3) a larger s and (s +s)/(ū +d) at small x (mainly due to LHCb W and Z rapidity data, and further enhanced by the ATLAS 7 TeV W/Z data in the CT18A and CT18Z fits). While the sensitivity of an individual tt data point can be similar to that of an individual jet data point at the LHC, the total sensitivity of the tt data is small due to the small number of tt data points. Hence, we did not find noticeable impact from the double differential distributions of the tt data included in the CT18 analysis. A similar finding was also reported in Ref. [102], in a CT14 analysis. To allow direct comparison to results obtained by the Lattice QCD (LQCD) community, we have also presented the CT18 predictions for PDF moments and sum rules in Sec. IV D, where we find good agreement between CT18 and LQCD results.
The final CT18 PDFs are presented in the form of 1 central and 58 Hessian eigenvector sets at NNLO and NLO. The 90% C.L. PDF uncertainties for physical observables can be estimated from these sets using the symmetric [23] or asymmetric [34,58]  For estimation of the combined PDF+α s uncertainty, we provide two additional best-fit sets for α s (M Z ) = 0.116 and 0.120. The 90% C.L. variation due to α s (M Z ) can be estimated as a one-half of the difference in predictions from the two α s sets. The PDF+α s uncertainty, at 90% C.L., and including correlations, can also be determined by adding the PDF uncertainty and α s uncertainty in quadrature [54]. Besides these general-purpose PDF sets, we provide a series of (N)NLO sets for α s (M Z ) = 0.111 − 0.123 and additional sets in heavy-quark schemes with up to 3, 4, and 6 active flavors.
Parametrizations for the CT18 PDF sets are distributed in a standalone form via the CTEQ-TEA website [168], or as a part of the LHAPDF6 library [169]. For backward compatibility with version 5.9.X of LHAPDF, our website also provides CT18 grids in the LHAPDF5 format, as well as an update for the CTEQ-TEA module of the LHAPDF5 library, which must be included during compilation to support calls of all eigenvector sets included with CT18 [170]. and (ii) the F p 2 , xF p 3 DIS structure function information extracted by CDHSW from ν-Fe data (Exp. IDs=108, 109).

ATLAS 7 TeV
Inclusive W/Z-production data. Regarding case (i) above, the ATL7ZW data set has been suggested as providing important information on the structure of the light-quark nucleon sea, and, in particular, favoring a larger-than-usual value for the strangeness suppression factor, R s (x, Q) of Eq. (10), i.e., larger than has typically been obtained in previous QCD analyses. Due to the dependence of the Drell-Yan cross section on antiquark PDFs, ∝q(x 1 )q(x 2 ), such data have often been suggested as playing an important role in disentangling the nucleon sea. Using the sensivity analysis and LM scans, we found that these data have a number of intriguing properties. While the correlation cosines of the ATL7ZW measurements with respect to various PDF flavors are quite modest (typically, | cos φ| < 0.6), the effect of including these data remains quite substantial, driving sizable shifts in the strangeness and other fitted PDFs, as we show in Sec. A 2 for CT18A NNLO.
The PDFSense and L 2 sensitivity methods discussed in Secs. II B 1 and V A 2 reveal an extremely pronounced pull of the ATL7ZW data on all light-antiquark flavors (s,ū, andd): see the L 2 -sensitivity profile for ATL7ZW in Fig. 61, as well in Fig. 73 for strangeness below.
This effect is largely attributable to the high experimental precision of the ATL7ZW data.
To examine how the inclusion of the ATL7ZW data into the original CT18 fit would modify the PDFs, we perform a global fit to produce the CT18A PDFs, in which the low luminosity (35 pb −1 ) sample of the same data (Exp. ID=268) was removed to avoid double counting. Our findings confirm those presented in Ref. [25], in which the ePump program was used to update the CT14 HERAII PDFs. The ATL7ZW data, with a total of 34 data points, cannot be fit well. Its χ 2 /N pt in the CT18A fit is found to be around 2.8, which is much larger than that for the full CT18A global data (about 1.191) with a total of 3674 data points. As a comparison, the full CT18 global fit yields χ 2 /N pt = 1.166 with a total of 3681 data points.
The ATL7ZW data have a sizable impact on the quark PDFs and their uncertainties. It is to decrease the u and d quark PDFs and increase the s quark PDF at x = 10 −4 ∼ 10 −3 , and increase the d PDF at x around 10 −2 and 0.3. Also, the error band of the d PDF is reduced significantly around x = 10 −3 , and the error band of the s-PDF is reduced for nearly all values of x. The large value of χ 2 /N pt associated with the ATL7ZW data indicates a significant discrepancy with respect to theory predictions that is not fully relieved by refitting. In Ref. [102], we investigated for possible tensions with other data sets included in the CT14 HERAII global fit (cf. Table I) by examining how the equivalent Gaussian variable S E of each data set varies when we increases the weight of ATL7ZW data in the ePump updating. We found that a few data sets were well-fitted before, but become poorly fitted after the weight of the ATL7ZW data is increased, the most significant ones being the NuTeVνµµ SIDIS (Exp. ID=125), the E866 σ pd /(2σ pp ) Drell-Yan data (Exp. ID=203) and the CMS 7 TeV electron asymmetry data (Exp. ID=267), cf. Fig. 25 of [102]. This can be explained by the sizable change found in the s-PDF andd/ū PDF ratio, cf. Fig. 26 of [102].
A similar conclusion also holds for the CT18A global fit. In the following sections, we will discuss more detailed comparisons of CT18A predictions and ATL7ZW data.
CDHSW data. In the final CT18Z fit, we decided to remove the CDHSW data (Exp. IDs=108, 109). It has been argued that these data sets were not analyzed properly and therefore should not be used in global fitting [171]. Furthermore, the LM scans, like the ones presented in Fig. 20, and the L 2 sensitivity plot in Fig. 24 indicate a preference of the CDHSW F 2 data (Exp. ID=108) for a harder g(x, Q = 100 GeV) at x ≈ 0.2, but a PDF Ratio to CT18NNLO x g(x,Q) at Q =100 GeV 90%C.L. CT18mCDHSW denotes the fit after removing the CDHSW data sets. The result of CT18Z fit is also shown for comparison.
softer gluon at x 0.5 -similarly to the CDF Run-2 jet data (Exp. ID=504). Therefore, one may wonder how the CT14 HERAII PDFs would change if we exclude these two CDHSW data sets from the original CT14 HERAII fit. Instead of redoing the whole global fit, we used ePump to quickly answer this question in Ref. [102]. There, we found that removing the CDHSW data led to slightly softer g-, s-,ū-, andd-PDFs at x > 0.1, and harder d-and u- PDFs at x > 0.03. Similarly, it was found that the Tevatron and LHC jet data also preferred a softer g-PDF at x > 0.2, cf. Fig. 6 of [102]. Thus, removing the CDHSW data was found to be more consistent with jet data.
In this work, we re-examine the effect of removing the two CDHSW data sets from a global fit. In this case, we remove the CDHSW data from the CT18 global analysis to a set of CT18-like PDFs, named "CT18mCDHSW" PDFs (i.e., CT18 "minus" CDHSW). In Given the pulls of the CDHSW data described above, a CT18-like global fit without these data can be expected to yield a softer gluon, relative to CT18, for x ∼ 0.2, but a hardening distribution at higher x > 0.5, just as is seen in the upper-left panel of Fig. 62. Excluding the CDHSW data has an opposing effect on the down-quark PDF in similar regions of x.

b. Alternative scale choices
Here, we examine the effect of changing the scales of DIS processes in the CT18 fit. This modified scale choice is the essential change explored in the CT18X alternate PDFs. As shown in Table III, the CT18X fit differs from the CT18 fit only by the change of the scales of DIS processes, which yields χ 2 /N pt = 1.148 with a total of 3681 data points.
Finally, we arrive at the CT18Z fit, in which we have made all those changes listed in  behavior which induces the slight reduction in the predicted Higgs production cross section from gluon-gluon fusion at the LHC as plotted in Fig. 50.
We observe a distinct pattern of PDF behaviors for the light-quark flavors under CT18Z/A/X, starting with the d-quark distributions shown in the panels of Fig. 64, which we present in an analogous fashion to the gluon PDFs of Fig. 63. In contrast to g(x, Q), the sensitivity of the Drell-Yan process to d(x, Q) leads to a moderate difference between the CT18 and CT18A NNLO results for d(x, Q), as seen in Fig. 64(a). This separation is realized as a mild, ∼ 1% suppression of the central CT18A distribution for d(x, Q) relative to CT18 about x ∼ 10 −3 , with a few-percent reduction in the accompanying PDF uncertainty.
The CT18Z fit in this region is pulled in the opposing direction, being enhanced by 2% with a comparable uncertainty for x < 0.  present in CT18X/Z, as can be deduced from the strongly similar PDF shapes obtained for these two fits in Fig. 63(b). Thus, the d-quark PDF (as well as the u) becomes softer at high x, while the gluon PDF becomes harder under CT18X/Z than with CT18 NNLO. The qualitative impact of fitting d(x, Q) at NLO, as opposed to NNLO, is fairly weak, being felt mostly at lower x 10 −3 ; again this effect is slightly amplified in the fits involving the modified scale choices, CT18Z/X NLO.
In Fig. 65, we examine the CT18Z/A/X NNLO alternative fits for theū-antiquark PDF, u(x, Q), in comparison with CT18 NNLO, plotting the PDF ratios with respect to the latter as before. Here, the remarkable property is the extent to which the deviations of CT18Z NNLO away from CT18 NNLO are driven by the modifications in CT18X, as attested by the very close agreement between the central PDFs and uncertainties for CT18Z/X shown in panel (b). This behavior in CT18Z recovered despite the weak suppression ofū(x, Q) in CT18A NNLO appearing in Fig. 65(a), which is almost entirely nulled once the DIS-and m c -scale choices found in CT18X are implemented on top of the ATL7ZW measurements.
We have noted already the important impact that the ATL7ZW data can have on the strange content of the proton, which remains poorly constrained compared to the light-quark PDFs. We illustrate the unique interplay of the CT18Z NNLO alternative fits for s(x, Q) in the panels Fig. 66, which follow the same conventions used in Fig. 65. The introduction of the ATL7ZW data leads to demonstrable deviations from the CT18 and CT18X NNLO fits This pattern of fitted strange distributions revealed in Fig. 66 can be interpreted in the context of the W , Z production cross section predictions shown in Fig. 51, in particular, the  In addition, we plot the strangeness suppression quantity, R s (x, Q), in Fig. 67. The qualitative behavior in this case is substantially driven by the patterns observed for s(x, Q) itself in Fig. 66. In addition, the rapid divergence at high x of the strangeness distribution is compensated in the flavor ratio.

120
W -boson and charm-jet (W + c) production from CMS at 7 TeV [113] in their fit, using an NLO calculation of MCFM [42]. The CMS 7 TeV W + c data are sensitive to s-PDF, which possibly helps relax tension with the ATL7ZW data. We do not include these CMS 7 TeV W + c measurements in CT18A(Z) because the full NNLO calculation is not yet available for them. In addition, the MMHT group [19] claims the tension between the ATL7ZW data with other sets may be relaxed by including NNLO calculations for DIS charm production as reported in Refs. [122,161]. In our own study, we confirm this trend, but find the improvement is too small to resolve the tensions we explore below. We also point out that MMHT adjusts the renormalization and factorization scales for their theoretical predictions of the ATLAS 7 TeV W/Z data, using µ R,F = M W,Z /2 [rather than µ R,F = M W,Z as chosen for CT18A(Z)], which may slightly affect χ 2 /N pt .
Theoretical predictions for ATL7ZW. It is useful to pursue the origin of these overall shifts in the description of the ATLAS W/Z data in the theory-to-data comparisons we plot in Fig. 68. As shown extensively in Sec. V B, we may explore the description of the ATL7ZW data under CT18A/Z as compared with CT18 NNLO. We do so in Figs. 68 and 69, plotting the rapidity distributions for W ± and Z production measured by ATLAS against our theoretical predictions in three panels of Fig. 68, and the corresponding values for the shifted residuals in Fig. 69, again for each channel. In Fig. 70  The tension in fitting the ATL7ZW data can be gauged by the ∆χ 2 -variations of the other CT18Z data sets in Fig. 71 (left). With the exception of the LHCb 8 TeV W/Z data (Exp. ID=250), the ∆χ 2 variations of all other plotted experiments are anti-correlated with ∆χ 2 for the ATLAS 7 TeV as the weight of these latter data is increased. As the ATLAS data are de-emphasized (w ∼ 0.1), the descriptions of the NuTeV SIDIS dimuonproduction data sets (Exp. IDs=124, 125) are most improved, further suggesting a strong tension between these measurements. This observation is consistent with the companion scan shown in Fig. 71 (right), which similarly plots the ∆χ 2 variations of the experiments most responsive to changing the weight of the NuTeV dimuon data. Especially considering N pt for the plotted experiments, the heavy over-weighing of the NuTeV data leads to a very rapid deterioration of χ 2 for the ATL7ZW points, which in fact worsens more quickly than the full CT18Z global analysis as the NuTeV weight is increased.
In addition, the growth in ∆χ 2 for other experiments as the weight of the ATL7ZW  in this data set. Among data sets with more modest numbers of points, the E866 ratio data (Exp. ID=203) show strong tension with ATL7ZW, in addition to the NuTeV dimuon data as well as the CMS 7 TeV µ-and e-asymmetry data (Exp. IDs=266, 267, respectively).
Experimental pulls on CT18Z PDFs. Regarding the ATL7ZW data, we find very strong pulls using the L 2 sensitivity approach upon the proton sea, i.e., the distributions of theū−,d−, and s-quarks. We have already noted in this Appendix the strong connection between the ATLAS data and the fitted strangeness distributions shown in Figs. 66 and 67.
We therefore highlight the competing pulls in CT18Z NNLO on R s (x, Q), which we show in and x = 0.1 (in the right). In both panels, we observe the strong preference of the ATL7ZW data for an enlarged R s , producing the enhancement seen in Fig. 67, especially for CT18A.
Also apparent in the panels of Fig. 72 are the opposing preferences of, e.g., the HERA and NuTeV SIDIS sets for a smaller value of R s , peaked more toward R s ∼ 0.5 at both values of x. We note that the pull of the ATL7ZW data is perhaps in greater conflict with the preferences of the other CT18Z data sets at higher x, all of which prefer a smaller R s than the full fit, with the exception of the F 2 CCFR neutrino data, which pull in the same direction as ATLAS, but more weakly. These findings are generally consistent with the L 2 sensitivity plot in Fig. 73, which shows the pulls of the CT18Z data on the strange-quark PDF, s(x, Q = 2 GeV).
Results for α s (M Z ) and m c . CT18Z involves a number of modifications that impinge upon inputs to perturbative QCD, especially the use of the x-dependent scale choice Q F,x .
This is due to the fact that the scale choice has a material effect on the description of the HERA I+II combined data, which in turn play a crucial role in determining the QCD evolution by merit of their broad coverage in Q. These differences are realized most clearly for m c increases slightly. In addition, the ATL7ZW data also exhibit a modest preference for larger charm masses, and these latter two experiments produce a small increase in the central preferred value of m c in CT18Z relative to CT18, although with a similarly significant uncertainty.
Appendix B: CT18 goodness-of-fit function and treatment of correlated uncertainties In this appendix, we summarize the implementation of the goodness-of-fit function χ 2 and marginalization of nuisance parameters in the CT18 family of fits. For the latter task, we normally follow a procedure, adopted since the CTEQ6 analysis, to estimate the correlated uncertainties using the correlation matrix published by the experiment. For the few data sets that do not provide the correlation matrix, we find it convenient to present the covariance matrix in an approximate form that separates the uncorrelated and correlated components.
The algorithm for this conversion is explained at the end of the appendix.
Definitions. In Eq. (1) of Sec. III A, we introduced the standard goodness-of-fit function χ 2 used in the recent CT fits. Here we review its treatment using a matrix notation.
We express D k , T k , and β kα in units of s k for each k. That is, we introduce a vector Similarly, λ ≡ {λ α } is a vector of length N λ ; and b ≡ S −1 β is a rectangular matrix of dimension N pt × N λ . In this matrix notation, Eq. (1) takes the form Minimization of χ 2 . The solution for the minimal value of χ 2 (with respect to nuisance parameters) is found in terms of the following matrices: Uppercase roman and script letters denote matrices of dimensions N pt × N pt and N λ × N λ , respectively. Therefore, I is an N pt × N pt identity matrix, and I is an For each theory input a, χ 2 can be minimized with respect to λ analytically, by solving for ∂χ 2 /∂λ α = 0 [23]. The solution is The global minimum χ 2 (a 0 , λ 0 ) for all experiments E can be found numerically as with d 0 ≡ S −1 (D − T (a 0 )), λ 0 ≡λ(a 0 ).
An equivalent form of this equation can be derived, where r 0 are the best-fit shifted residuals: The representation (B8) is particularly informative. We anticipate that, in a good fit of theory to an experiment E, the shifted residuals r 0 , quantifying agreement with individual data points, as well as the nuisance parametersλ 0 , quantifying the systematic shifts, are distributed according to their own standard normal distributions, N (0, 1). Comparisons of the two empirical distributions to the expected N (0, 1) distributions serve as the tests for the goodness of fit and for the implementation of systematic errors [33]. 128 In particular, a diagonal element (cov) ii [no summation] is the quadrature sum of the statistical, uncorrelated systematic, and correlated systematic uncertainties: (cov) ii = s 2 i,stat + s 2 i,uncor sys + s 2 i,cor sys , where s 2 i,cor sys ≡ α β 2 iα . With the help of Eq. (B10), the shifted residuals in Eq. (B9) then are written as Computing them thus requires that we know the full uncorrelated error s i = s 2 i,stat + s 2 i,uncor sys . Finding a correlation matrix from the covariance matrix. In some experimental measurements, such as the LHCb 8 TeV W/Z production [94], only the full covariance matrix is provided, making impossible the straightforward reconstruction of the shifted residuals according to Eq. (B14). In several cases, we find it feasible to iteratively reconstruct the approximate uncorrelated systematic and correlated systematic contributions, s k,uncor sys and β kα , from the provided covariance matrix cov ≡ K 0 , by assuming that the systematic shifts are dominated by a certain number M λ , with M λ ≤ N λ , of fully correlated linear combinations. The approximation makes use of the positive-definiteness of the covariance matrix and its diagonal elements, cf. Eq. (B13). It represents the original covariance matrix by a numerically close matrix given by the sum of a diagonal matrix Σ, interpreted as consisting of total uncorrelated errors, and a non-diagonal square matrix K, interpreted as a product of the correlation matrix β and its transpose.
In particular, suppose we find the eigenvalues x 2 k of K 0 and sort them in the descending order: Here O is an orthogonal matrix. We partition x into a matrix y containing the largest M λ eigenvalues and a matrix z with the smallest (N pt − M λ ) ones: y ≡ diag x 2 1 , x 2 2 , ..., x 2 M λ , 0, ..., 0 , z ≡ diag 0, ..., x 2 M λ+1 , ..., x 2 Npt .
Recalling that the diagonal elements of matrices Y = O T yO and Z = O T zO are nonnegative, we then construct a diagonal matrix Σ 1 ≡ diag Z 11 , ..., Z NptNpt with Z ii > 0, and another positive-definite matrix, K 1 ≡ Y + Z − Σ 1 . We can iterate the steps in Eqs. (B15)-(B18) by computing Σ a+1 and K a+1 at step a as Σ a+1 ≡ Σ a + diag Z 11 , ..., Z NptNpt , (B19) Here the matrices Y and Z are recomputed in each step using K a as the input. After a sufficient number of steps a, the sum approaches an asymptotic matrix that is close to the input matrix K 0 in the sense of the L p norm Npt i,j=1 |(K 0 ) ij − (C a ) ij | p with p = 2 or 1. If the extraction of the uncorrelated and fully correlated components is feasible, the asymptotic L p distance can be made small by choosing a large enough M λ . For the three experiments [91,93,94] in the CT18 data set that provide only the covariance matrices, the M λ values giving good convergence lie in the range between N λ /2 and N λ .
By comparing Eqs. (B11) and (B21), we identify Hence, s i is estimated from (Σ a ) ij ; and β iα can be estimated from (K a ) ij by singular values decomposition.

130
Appendix C: Non-perturbative parametrization forms As noted in Sec. III C, the philosophy of the CT global analyses is to explore a broad range of parametric forms for the starting-scale distributions at Q = Q 0 . The goals of these investigations are (i) to select a sufficiently flexible functional form capable of fitting an expansive high-energy data set without overfitting; and (ii) to understand the uncertainties associated with the choice of parametrization.
In CT18, we use a nonperturbative parametrization generally similar to the one that served as the basis for CT14 HERAII , but with a slightly more flexible parametrization for the sea-quark distributions. This enhanced flexibility is necessitated by the inclusion of LHC Run-1 data with direct sensitivity to theū,d and strangeness content of the nucleon. The functional form used to parametrize the starting-scale valence distributions of the u-and d-quark (i.e., u v and d v ) is a polynomial in the parameter y ≡ √ x, P v a (y) = a 3 (1 − y) 4 + a 4 4y(1 − y) 3 + a 5 6y 2 (1 − y) 2 + a 6 4y 3 (1 − y) + y 4 , where a 6 = 1 + 1 2 a 1 .
Owing to the comparative lack of empirical constraints to nucleon strangeness, not all the parameters above are permitted to float freely for s(x, Q) in CT18(Z); rather, we constrain a 4 = a 5 and a 6 = a 7 for the strangeness parameters. As with the valence distributions, we constrain the prefactor exponents as to ensure a finite value for the strangeness suppression ratio, R s = [s +s] [ū +d], in the limit x → 0, and the convergence of . We note that the final normalizations of the PDFs are determined by the parametric forms reported here as well as additional overall factors implemented to guarantee the momentum and flavor-number sum rules. The ratio of these overall normalization factors governs the specific finite value of R s in the x → 0 limit.
In Table X, we summarize the central fitted values of the parameters noted above for CT18 (upper rows) and CT18Z (lower rows). Those parameters entries marked with "−" do not participate as degrees of freedom for the associated flavor.  The functional forms associated with each parametrization are defined explicitly in the text of this section, Eqs. (C1)-(C7). Those entries corresponding to parameters that are not actively fitted for a given PDF flavor are indicated by a dash, "-." For the sea-quark distributions,ū,d, and s, we fix a 3 = 4, and for s alone, we fix a 8 = 1, as given parenthetically above for CT18(Z).

Appendix D: Fitting code developments
The inclusion of LHC data requires further improvement on efficiency of the global analysis machinery. The usage of various fast interface on the calculations of structure functions and cross sections becomes mandatory and conventional. However, the global analysis remains a time consuming task due to the large numbers of both experimental data and input parameters. Figure 76 shows the full structure of a typical global analysis taking CT18 as an example. The structure can be separated into three major layers. LY0 represents the out layer of the analysis that repeats the fit with different inputs or constraints, e.g., QCD coupling or heavy quark masses for studying parametric uncertainties, non-perturbative parametrization forms for parametrization dependence, and Lagrange Multiplier scans on certain observables or scans on weight of certain data set. The middle layer, LY1 of the actual fit, scans over the PDF parameter space and constructs its χ 2 profile. That provides the best-fit PDF and the associated uncertainties either in a form of Hessian PDFs as in CT18 or MC replicas. The calculation of global χ 2 , LY2, is the core part of the fit that will be repeated for every point in the scan of the PDF parameters. It consists of calculations of the cross sections for thousands of data points and computation of individual χ 2 for each data set.
Most computational efforts shown in Figure 76 can be easily parallelized. For instance, in LY0, one can simply submit thousands of variation fits at the same time to a large computing cluster since those fits are mostly independent. In LY1, the major task is to find a global minimum of the χ 2 in a parameter space with large dimensions (N para. ∼ 30). The possible parallelization depends strongly on the minimization algorithm. In CT18 analysis we use the variable metric method as implemented in the MINUIT package [172]. It involves numerical calculations of the first-order derivatives of the χ 2 , second-order derivatives, and sequential minimum searches along fixed directions, in the PDF parameter space. The former two calculations are highly parallelizable with a scaling of 1/N para. and 2/(N 2 para. + N para. ) respectively. For the core part, again the calculation of individual χ 2 of different data sets, including their cross sections, can be done simultaneously.
Concerning the realization at the code-programming level, the solution can be either using MPI or OpenMP. The latter are restricted to platforms with shared memory, but requires much less efforts on adjusting our fitting code that was previously intended for sequential run only. We use a similar approach as OpenMP that minimizes the changes inside our fitting code. Being specific, we use the fork Linux command to branch the main program into multi-threads, each with a copy of the master memory and carrying out the calculations independently. Later we use the join command to collect and return the individual results to the main program. The implementation of this fork-join algorithm is borrowed from the widely-used CUBA library [173] for multi-thread Monte-Carlo integration.  computed here using the CT14 HERAII fit under several error treatment scenarios. χ 2 /N pt is first given without implementing any decorrelation scheme ("original data"); using the experimentally-proposed decorrelation scheme described above ("+ decorr."); and adding on top of these decorrelations an overall uncorrelated Monte Carlo uncertainty of 0.5% in the final column. and otherwise L = 0 for x < x min or L = 1 for x > x max . We note that Eq. (E1) governs the magnitudes of the decorrelated uncertainties, which otherwise inherit the sign of the original uncertainty. A similar algorithm is applied to decorrelate the JES flavor response JES62.
The result of applying these options is a 90 unit reduction in the total χ 2 of the ATLAS jet data consistent with the summary given in Table XI below. In addition, we have explored the effect of including an additional overall uncorrelated uncertainty ('MC unc.') to accommodate the estimated intrinsic statistical noise associated with Monte Carlo generation of NNLO/NLO K-factors. The incorporation of these additional uncertainties leads to the further reductions in χ 2 shown in Table XI.
The xFitter-profiled PDFs also include second-order diagonal terms, ∼ 1 2 (λ min k,th ) 2 (f + k + f − k − 2f 0 ). We note that the off-diagonal terms, ∼ λ min k,th λ min j,th for k = j, are not presently included since the off-diagonal, second-order partial derivatives cannot be constructed solely in terms of Hessian PDF error sets [139].
Here, we use the Hessian-profiling method of xFitter, as well as ePump updating, to explore the impact of the ATLAS 7 TeV inclusive W/Z-production data [35] on several PDF sets. The change in the total χ 2 values before and after profiling/updating with the ATLAS 7 TeV Z/W -production data is presented in Table XII for each PDF set, calculated within both the xFitter and ePump programs. First, we explore all 7 measurements, having a combined N pt = 61 data points: W + , W − , neutral current DY in the low-mass, Z-peak, and high-mass regions for the central and forward selections. In this case, we can directly compare with the ATLAS [35] and MMHT19 [19] analyses. As a second case, we take only the 3 most precise measurements, i.e., the W + , W − , and Z-peak DY data for the central selection, with N pt = 34 data points in total, which are included in the NNPDF3.1 [3] and CT18A(Z) global analyses. The comparison of the fitted χ 2 /N pt values for these data in CT18A(Z), NNPDF3.1 and MMHT19 can be found in Sec. A 3.  well with the results presented in Ref. [35]. Here, we should apply the tolerance T = 1.645 in xFitter as the CT PDFs are defined according to a 90% C.L. [1,29]. As shown in Ref. [25], the same results can be reproduced by the ePump code when the tolerance is set to T = 1.645.
However, as clearly discussed in [25], setting T = 1.645 in the ePump calculation for the CT PDFs is equivalent to assigning a very large weight (about 100/1.645 2 ) to the new data set included in the fit. xFitter profiling therefore generally overestimates the impact of new data sets when using the CT PDFs. Meanwhile, a universal tolerance value (T = 1.645) is not able to capture the constraints of Tier-2 penalty to determine the CT PDF error sets. An appropriate way to update an existing CT PDF set with the inclusion of any new experimental data is to adopt a dynamical tolerance. As such, xFitter profiling yields a smaller χ 2 value (63) than does ePump updating with a dynamical tolerance (144), as can be seen by comparing the rightmost entries in the last row of Table XII. This conclusion also holds when using xFitter profiling with the MMHT [2] and PDF4LHC15 [52] PDFs.
The χ 2 value of the 34 highest-precision ATLAS 7 TeV W/Z-production data points is found to be χ 2 = 87.6 in the CT18A global fit, cf. Table IX. To this we compare the corresponding value found using ePump updating with dynamical tolerance, for which we obtain χ 2 = 144, as reported in Table XII. This value includes two district contributions, from the difference between theory and data of the ATLAS 7 TeV W/Z-production itself as well as from theoretical nuisance parameters, which can be interpreted as the increase of the "old" data set included in the CT18 fit. We found that these are 104 and 40, respectively.
The large increase in χ 2 originated from the "old" data set, after ePump updating, indicates the presence of some tensions between the ATLAS 7 TeV W/Z-production data and the other data sets included in the original CT18 fit, such as the SIDIS dimuon data [25], cf. Fig. 71.
These differences between the χ 2 values in CT18A (87.6) and from ePump updating with a dynamical tolerance (104) indicates the breakdown of the linear approximation used in the Hessian-updating method, when applied to this case. This can be confirmed by examining the updated central PDF set. In Fig. 77, we show the gluon and strangeness PDFs at Q = 100 GeV for the CT18 (with 90% C.L. error) and CT18A global fits, compared with ePump dynamical-tolerance updating and xFitter profiling of CT18, based on the inclusion of the ATLAS 7 TeV W/Z-production data. As compared to the CT18 PDFs, we see that the updated g and s PDFs from the ePump program are similar to the CT18A PDFs, although with somewhat smaller shifts in the data-sensitive range, 10 −3 < x < 10 −1 . In contrast, the default xFitter profiling produces a much larger shift in both g-and s-PDFs, than the CT18A global fit. As a result, xFitter profiling produces a too large change in the s-PDF so that its central prediction touches the upper error band of CT18 at x > 0.02. We have found similar features in the comparison of other flavor PDFs. Namely, the default xFitter program generally overestimates the impact of new data sets in updating the existing PDFs [25].

Comparison of different NNLO predictions.
In addition to the studies described above, we have also used the xFitter and ePump frameworks to examine aspects of the theory calculations for the ATLAS 7 TeV W/Z data.
In the CT18A(Z) global fits, the NNLO predictions for these measurements were calculated using NNLO/NLO K-factors implemented via APPLgrid files. Specifically, the K-factors used in our CT18A(Z) fits were directly extracted from xFitter, which were calculated with the DYNNLO code [40,41]. It was noted by the ATLAS Collaboration in Ref. [35] that the integrated fiducial Z , W + and W − cross sections predicted by the NNLO codes FEWZ [45][46][47] and DYNNLO [40,41] differ by about 0.2%, 1.2% and 0.7%, respectively. We have directly confirmed this result in our own calculations. measurements of the ATLAS 7 TeV W/Z-production data. First, we verified that the NLO predictions agree, within 0.2%, among DYNNLO, MCFM [43,44] and FEWZ, which is the result of the fact that all the three codes take the dipole formalism [177,178] in the NLO calculations. Since the NNLO predictions can be expressed as the NLO predictions multiplied by their K-factors, we compare in Fig. 78 the K-factors obtained from each code, for every ATLAS 7 TeV W/Z-production data point (with 34 data points in total). It can be seen that the differences among these three K-factors varies as a function of the Z-boson rapidity, or, the rapidity of the charged lepton from the W -boson decay. These differences can be sizable, 1%, as compared to the typically sub-percent statistical uncertainty found in the ATLAS 7 TeV W/Z data. We also find that the predictions of MCFM generally lies between those of DYNNLO and FEWZ. This difference can be understood as a consequence of the the different NNLO schemes: FEWZ adopts sector decomposition [179][180][181], while DYNNLO and MCFM take transverse-momentum (q T ) [40,41] and N -jettiness (T N ) [43] subtractions, respectively.
It is also useful to investigate the dependence of the ATLAS 7 TeV W/Z fit quality on the specific choice of NNLO calculation scheme. We have carried out such a study, and summarize our findings in