Impact of LHC vector boson production in heavy ion collisions on strange PDFs

The extraction of the strange quark parton distribution function (PDF) poses a long-standing puzzle. Measurements from neutrino-nucleus deep inelastic scattering (DIS) experiments suggest the strange quark is suppressed compared to the light sea quarks, while recent studies of W/Z boson production at the LHC imply a larger strange component at small x values. As the parton flavor determination in the proton depends on nuclear corrections, e.g. from heavy-target DIS, LHC heavy ion measurements can provide a distinct perspective to help clarify this situation. In this investigation we extend the nCTEQ15 nPDFs to study the impact of the LHC proton-lead W/Z production data on both the flavor differentiation and nuclear corrections. This complementary data set provides new insights on both the LHC W/Z proton analyses and the neutrino-nucleus DIS data. We identify these new nPDFs as nCTEQ15WZ. Our calculations are performed using a new implementation of the nCTEQ code (nCTEQ++) based on C++ which enables us to easily interface to external programs such as HOPPET, APPLgrid and MCFM. Our results indicate that, as suggested by the proton data, the small x nuclear strange sea appears larger than previously expected, even when the normalization of the W/Z data is accommodated in the fit. Extending the nCTEQ15 analysis to include LHC W/Z data represents an important step as we advance toward the next generation of nPDFs.


I. INTRODUCTION
Parton distribution functions (PDFs) are key elements required to generate concrete predictions for processes with hadronic initial states in the context of QCD factorization theorems. The success of this theoretical framework has been extensively demonstrated in fixedtarget and collider experiments (e.g., at the TeVatron, SLAC, HERA, RHIC, LHC), and will be essential for making predictions for future facilities (EIC, LHeC, FCC). Despite the above achievements, there is yet much to learn about the hadronic structure and the detailed composition of the PDFs [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18]. Contribution of strange initiated channels to W ± and Z boson production for proton-lead (pPb) at the LHC. The blue lines represent total cross-sections, the yellow lines are cross-sections with the strange initiated channels subtracted. The lower panels show the ratio compared to the total cross-section.
Although the up and down PDF flavors are generally well-determined across much of the partonic x range, there is significant uncertainty in the strange component, s(x). The strange PDF is especially challenging because, in many processes, it is difficult to separate it from the larger down component. However, as we push to higher precision and energies, an accurate determination of the strange PDF is, apart from its intrinsic fundamental importance, essential not only for LHC measurements, but for a wide variety of processes [1,[11][12][13][14][15][16][17][18]. For example, the knowledge of the nuclear strange distribution in heavy nuclei is crucial for providing a reliable baseline for hard probes of the quark gluon plasma (QGP) which is characterized by enhanced production of strangeness [19][20][21]. Additionally, small x nuclear PDFs are essential for computing the composition of air showers from ultra-high energy cosmic rays [22][23][24][25].
The recent results from the LHC for W/Z boson production in pp collisions predict a large strange to lightsea ratio [26][27][28][29][30][31]. This is a rather surprising result as it differs from earlier determinations based on analyses of neutrino deep inelastic scattering (DIS) data from NuTeV and CCFR experiments [32][33][34] or charged kaon production data from HERMES [35]. See Ref. [36] for more details on the earlier determinations of the strange distribution, and Ref. [37] for a study of the compatibility of the ATLAS and CMS results using the xFitter framework [38].
It is not easy to directly compare the pp LHC results with the fixed target experiments, as the earlier measurements were generally done using nuclear targets (typically Fe or Pb). Additional complications arise from the fact that there is a controversy about the proper nuclear correction factors for the charged current (CC) and neutral current (NC) DIS measurements [39][40][41][42][43]. As a result, the choice of heavy target neutrino DIS data sets varies widely among not just the many nuclear PDF (nPDF) determinations, but also for the proton PDF fits [1,2]. Moreover, in the proton case the nuclear corrections are applied in different ways.
Conversely, it was already demonstrated that the W/Z LHC data can provide some important information on the strange and gluon nPDFs [4,6,44]. To demonstrate the impact of the heavy ion W ± /Z data on the strange PDF, in Fig. 1 we display the contribution of the strangeinitiated process as a function of rapidity. We observe the strange component can be as much as 20% to 30% of the total. For this reason, we concentrate in the following on the constraints for the nuclear strange and gluon distributions given by the W and Z data from proton-lead (pPb) collisions at the LHC. This process is an ideal QCD "laboratory" as it is sensitive to i) the heavy flavor components {s, c, ...}, ii) the nuclear corrections, and iii) the underlying "base" proton PDFs. Such an analysis provides an independent perspective on the subject and can help disentangle the flavor separation and nuclear modifications.
In the current investigation, we will study the production of W and Z bosons in proton-lead (pPb) collisions at the LHC; this involves similar considerations as the pp case, but also brings in the nuclear corrections. We will be focusing, in particular, on the strange and gluon distributions to see how these are modified when the LHC measurements are included. In Sec. II we review the various data sets used in our analysis along with the separate fits extracted. In Sec. III we present the quality of the fits and comparisons of data with the theory, and demonstrate the impact on the resulting PDFs. In Sec. IV we compare our final PDF fit with other results from the literature. In Sec. V we recap the key outcomes of this study.

A. The nCTEQ++ Framework
The nCTEQ project extends the proton PDF global fitting effort by fully including the nuclear dimension. 1 Previous to the nCTEQ effort, nuclear data was "corrected" to isoscalar data and added to the proton PDF fit without any uncertainties [45]. In contrast, the nCTEQ framework allows full communication between the nuclear data and the proton data; this enables us to investigate if observed tensions between data sets could potentially be attributed to the nuclear corrections.
The details of the nCTEQ15 nPDFs are presented in Ref. [3]. The present analysis is performed in a new C++ code base (nCTEQ++) which enabled us to easily interface to external programs such as HOPPET [46], APPLgrid [47], and MCFM [48]. The nCTEQ15 fit has been reproduced in this new nCTEQ++ framework.
For the current set of fits, we use the same 16 free parameters as for the nCTEQ15 set, and additionally open up three parameters for the strange PDF, for a total of 19 parameters. Recall that for the nCTEQ15 set, the strange PDF was constrained by the relation s =s = (κ/2)(ū +d) at the initial scale Q 0 = 1.3 GeV so that it had the same form as the sea quarks.
Our PDFs are parameterized as where k = {1, ..., 5}. The 16 free parameters used for the nCTEQ15 set model the x-dependence of the {g, u v , d v ,d +ū} PDF combinations, and we do not vary thed/ū parameters; see Ref. [3] for details. To this, we now add three strange PDF parameters: {c s+s 0,1 , c s+s 1,1 , c s+s 2,1 }; these parameters describe, correspondingly, the overall normalization, the low-x exponent and the large x exponent of the strange distribution.

B. Experimental Data Sets
In this analysis we use the deep inelastic scattering (DIS), Drell-Yan (DY) lepton pair production, and RHIC pion data employed in our earlier nCTEQ15 analysis [3]. Additionally, we use W and Z inclusive data from protonlead collisions at the LHC. Specifically, we include the following data sets: ALICE W ± boson production [54,55], ATLAS Z boson production [50], ATLAS W ± boson production [49], CMS Z boson production [52], CMS W ± boson production [51], CMS Run II W ± boson production [53], and LHCb Z boson production [56]. The data sets are outlined in Table I

TABLE III
We present the χ 2 /N dof for the individual data sets, the individual processes {DIS, DY, Pion, LHC}, and the total. We also show the total χ 2 contribution from the LHC normalization penalty. Note that the Pion χ 2 /N dof is shown for comparison, but this data is only included in the nCTEQ15WZ fit.
calculations are performed at the next-to-leading order (NLO) of QCD. In particular, the calculations for W and Z boson production have been performed using the MCFM-6.8 program [57] interfaced with APPLgrid [47] in order to speed up the computations.

C. The PDF Fits
We now use the nCTEQ++ framework to include the LHC W ± /Z pPb data and extend the nCTEQ15 fit. Comparing the LHC pPb data to the nCTEQ15 results, we find that these data generally lie above the theory predictions [44]; hence, if we allow for a normalization uncertainty, this additional freedom can significantly improve the fit. The experiments have an associated luminosity uncertainty (cf., Table II), and we will use this as a gauge as we shift the data normalizations. It is reasonable to tie the normalizations for {W ± , Z} data from individual experiments (e.g, CMS Run I) to a single normalization factor as these uncertainties are fully correlated.
Previous studies implied a close connection between the normalization of the W ± /Z data and the extracted strange PDF [44]. To systematically investigate the effect of the normalization in detail, we will use a series of fits outlined in Table III and summarized below: nCTEQ15: This is the original set of nuclear PDFs as computed in Ref. [3].
Norm0: We include the LHC pPb data, but we do not allow for any floating normalization of the LHC data.
Norm2: We include the LHC pPb data, and allow for 2 normalization factors; one for the ATLAS Run I data, one for CMS Run I; we do not renormalize CMS Run II data in this fit.
Norm3: We include the LHC pPb data, and we allow for 3 normalization factors; one for the ATLAS Run I data, one for the CMS Run I data, and a separate one for the CMS Run II data.
nCTEQ15WZ: This is the same as Norm3, but we also include the RHIC inclusive pion data directly in the fit. This is discussed in Sec. IV.
All four of these new PDF fits are based on the DIS and DY data from the nCTEQ15 analysis and the LHC data sets, as outlined in Sec. II B and Table I. As with our nCTEQ15 study, we will present results both with and without the inclusive pion data [58,59]. For the comparison of the the W ± /Z normalizations fits {Norm0, Norm2, Norm3}, we will not include the pion data; however, we do compute the pion χ 2 , as shown in Table III, to demonstrate the compatibility. 2 In Sec. IV we then present a separate fit, nCTEQ15WZ, which does include the pion data. As we will see the impact of the pion data is marginal.
Normalization Factors: Table II shows the determined normalization factors used in each fit. All the normalization shifts are between 1σ and 2σ of the quoted normalization uncertainty, but all are systematically below unity; the appropriate normalization penalties are included in the χ 2 calculations. The detailed prescription we use for fitting data normalizations is provided in the Appendix A. For the ALICE and LHCb sets, the current uncertainties provide sufficient flexibility that we do not use an additional normalization factor for these data.
Quality of individual data sets: To provide more details regarding the source of the χ 2 contributions, in The χ 2 /N dof bar charts provide incisive information as to which data sets are driving the fit. We discuss each fit in turn.
nCTEQ15: Starting with the nCTEQ15 set, we note that (except for a few outliers) the DIS and DY data is well described by these PDFs; 5 by comparison, the LHC W ± /Z data (which was not included in the original nCTEQ15 fit) is not well described. As was detailed in Ref. [44], an important contribution to this large χ 2 comes from the small x region where the nuclear PDFs are poorly constrained. The re-weighting analysis of Ref. [44] demonstrated that we can improve the fit by adjusting the small x behavior of the PDFs, but this alone will not bring all the data sets into the range of χ 2 /N dof ∼ 1; something else is required.
Norm0: As a first step, this fit includes the LHC W ± /Z data, but does not include any floating normalization factors. This fit will tell us the extent to which we can adjust the PDFs to fit the LHC data before we begin to adjust the normalization factors. Examining Fig. 2, we see the impact of this fit on the DIS and DY data is generally small for many of the data sets, but does result in noticeable improvement for a few of the sets including 5115 (NMC Ca/D) and 5121 (NMC Li/D), for example. 3 When we refer to N dof for the total χ 2 we calculate it in the usual way as the difference between number of data points and number of free parameters (N dof = N data − Npar). However, when referring to N dof for individual experiments or data sets we set it to be equal to the number of data points (N dof = N data ). 4 The IDs of the specific non-LHC experiments can be found in Ref. [3]. In general, DIS data sets are 51xx, DY sets are 52xx, and W ± /Z sets are 62xx. 5 We find the DIS experiments 5108 (Sn/D EMC-1988) and 5121 (Ca/D NMC-1995) to be outliers with χ 2 /N dof > 2; this is consistent with other analyses [4,60].
However, it does significantly improve the LHC W ± /Z fit reducing the partial χ 2 /N dof of this data from 6.20 to 1.47 for the 120 LHC data points. Although this is a notable improvement, a number of the LHC data sets still have χ 2 /N dof values well above one.
Norm2: In this fit, we allow two floating normalization factors (one for ATLAS and one for CMS Run I) which are allowed to vary in the fit. The contribution of the normalization penalty is included in the total χ 2 . We see the impact of the floating normalization factors on the DIS and DY data is again small, as was the case for the Norm0 fit. But the Norm2 fit dramatically improves the LHC W ± /Z data reducing χ 2 /N dof of this data to 1.15 as compared to 1.47 for the Norm0 fit. While the Norm2 fit is a substantial improvement over the Norm0 results and all LHC data sets have χ 2 /N dof < 2, there are still a few sets at the upper limit of this range.
Norm3: Finally, we now perform a fit with three normalization factors: one for ATLAS (Run I), and one for each CMS Run I and CMS Run II.
As before, the modifications to the DIS and DY sets are minimal, but we do continue to see an improvement in the LHC sets; namely the χ 2 /N dof of this data improves to 0.91 as compared to 1.15 for the Norm2 fit.
Comparing Norm3 with nCTEQ15 for the other data sets, we see that the χ 2 /N dof for the DIS data is essentially the same (0.91 vs. 0.90), the DY increases slightly (0.73 vs. 0.77), and the pion χ 2 (computed a posteriori) increases slightly as well (0.25 vs. 0.39); these differences are relatively small compared to the significant improvement in the LHC data (6.20 vs. 0.90).

B. Comparison of Data with Theory
To obtain a more complete view of the fit quality, in Figs. 3, 4, 5, and 6 we display the comparison of the LHC data with theory predictions. The data points and errors are taken directly from the experimental measurements. However, it is important to note that we have shifted the theoretical predictions by the appropriate normalization factors; this allows us to present the fits with different normalizations on a single plot, and provides a more accurate visual description of the quality of the fit.

Large x Region:
Our first observation is that the experimental data consistently lies above our theoretical predictions. From Table II, recall that all the fitted normalization factors are less than one, indicating that the fit prefers a reduction of the data values, typically in the range of ∼ 5%; because we have shifted the theory, this is not as obvious in Figs. 3-6.    Even with the normalization shifts, we see the theory predictions still lie well below the data for a number of sets. This is most evident in the negative y region for the Run I W − data sets 6211 (ATLAS W − ) and 6231 (CMS W − Run I), and to a lesser extent 6215 (ATLAS Z). Interestingly, the Run II data generally show good agreement across the full y range.
The negative rapidity region corresponds to the large x region of the lead PDF. The large x region is already rather well constrained by the fixed-target measurements, so there are limits as to how much the new LHC data can shift the PDFs in this region. Also note, that in the large x region we are in the "anti-shadowing region" (x ∼ 0.1) where the nuclear corrections typically enhance the nuclear PDF relative to the proton. Thus, not including the nuclear corrections in this region would increase the discrepancy.

Small x Region:
In the large rapidity (small x) region, we generally find good agreement between our new fits and the data. But, this is in striking contrast to the nCTEQ15 PDF which lies well below many of the data points at large y; this behavior is clearly evident, for example, in 6215 (ATLAS Z), and 6213, 6232, 6233, 6234 (CMS W ± Run I and Run II). Clearly, the new LHC W ± /Z data provides important new PDF constraints in this kinematic region that were not available in the nCTEQ15 analysis.
As larger rapidity corresponds to smaller x values, this puts us in the "shadowing region" (x < ∼ 0.1) where the nuclear PDFs are generally expected to be suppressed relative to the proton. If the nuclear shadowing correction were reduced in this region, that would bring the theory closer in line with the data without the need for large normalization factors. The precise value of the nuclear corrections is still an open question; for example, Refs. [39,40,61] found that the shadowing correction for the νN charged-current neutrino DIS was reduced as compared to the ± N neutral-current DIS. If such an adjustment were applied to the LHC W ± /Z data, it would move the theory closer to the data and reduce the normalization factor. Disentangling the nuclear effects from the underlying parton flavor components is intricate, and a reanalysis of the neutrino DIS data is currently in progress [62].

C. The PDFs
Finally, we make a detailed examination of the underlying flavor PDFs from these various fits. In Figs. 7-9 we display the nPDFs for a full lead nucleus at three separate scales. The lowest scale (Q = 2 GeV) is close to our initial evolution scale of Q 0 = 1.3 GeV, the largest scale (Q = 90 GeV) is in the range relevant for W ± /Z production, and the intermediate scale (Q = 10 GeV) helps illustrate the effects of the DGLAP evolution.
We choose to display the full lead nPDF as this is the physical quantity which enters the calculation. 6 This is computed using: and we assume isospin symmetry to derive the neutron PDF.

Strange and Gluon nPDFs:
Examining the curves for up and down distributions, we see there is minimal variation between different fits as these flavors are strongly constrained by other data. Interestingly, we also see that the small x uncertainty is reduced at higher scales (see Figs. 8 and 9). We observe a slight modification in theū andd distributions as these are closely linked to the gluon and strange distributions which we will discuss in the following.
Turning to the gluon and strange PDFs, we see significant differences. In particular, the fits seems to prefer a larger value for both the gluon and strange PDFs at intermediate x values, which is the region relevant for the LHC heavy ion W ± /Z production. We discuss these fits in turn.
Norm0: Examining the Norm0 fit for Q = 2 GeV (Fig. 7), we see a distinct excess in the strange and gluon PDFs in the region x ∼ 0.03; this is also evident in Fig. 10 where we have plotted the ratio relative to the nCTEQ15 values. At Q = 2 GeV, the peak of the gluon and strange  Recall that the Norm0 fit does not allow any normalization adjustment in the fit. Since the data consistently lie above the theoretical predictions, it appears that the Norm0 fit is exploiting the uncertainty of gluon and strange PDFs to try and pull up the theoretical predictions in line with the data by increasing the PDFs in the relevant x region. Additionally, we observe a similar (but less pronounced) behavior in theū andd distributions. As momentum must be conserved, we see the Norm0 strange PDF dips below nCTEQ15 at both high and low x values, while the gluon is below nCTEQ15 at higher x values. Part of the reason the deformation of the gluon and strange PDFs is so large at Q = 2 GeV is to compensate for the DGLAP evolution which will tend to diffuse the excess in the gluon and strange distributions at the Q = 90 GeV scale, cf., Fig. 9.
Norm2 and Norm3: In contrast to the Norm0 result above, the Norm2 and Norm3 fits allow us to investigate the effect of including the normalization parameters into the fit; this is crucial in reducing the χ 2 /N dof for the LHC heavy ion data. The effect on the resulting nPDFs is evident as shown in Fig. 7 where we see that the excess in both the strange and gluon is systematically reduced as we introduce normalization parameters. In Fig. 7 we also observe the greatly increased error band on the Norm3 strange PDF as compared to nCTEQ15; this is of course due to the additional fitting parameters for the strange quark included in the Norm3 analysis.
To highlight the magnitude of these differences, in Fig. 10 we plot the ratios of the PDFs compared to nCTEQ15. At Q = 2 GeV, we see that the Norm0 gluon is nearly a factor of 2 times the nCTEQ15 value, with a peak at x ∼ 0.03. The Norm2 and Norm3 gluon PDFs are reduced to ∼ 60% and ∼ 40% above nCTEQ15, respectively. 8 Similarly, at x ∼ 0.03 and Q = 2 GeV, the strange PDF for both Norm0 and Norm2 are ∼ 60% above the nCTEQ15 value, while the Norm3 result is reduced to ∼ 25%. In Fig. 10 we also display ratio at Q = 90 GeV which illustrates the effect of the DGLAP evolution. We see that the gluon is now reduced to ∼ 15% above the nCTEQ15 value, the strange is reduced to ∼ 25% above the nCTEQ15 value, and both peaks have shifted to lower x values.
Because the DGLAP evolution has "washed out" the detailed peak structure at low Q values, it is necessary for the fit to amplify the distortion at low Q so that a remnant of the effect survives at high Q. Nevertheless, the remaining excess at Q = 90 GeV is sufficient to improve the χ 2 of the fits.
Additionally, we note that the heavy-flavor reweighting analysis of Ref. [63] also observed an increase of the gluon 8 Note we are focusing here on the intermediate x region (x > ∼ 0.01) not only because this is the central x region for W ± /Z production, but because the small x region is poorly constrained. nPDF in the intermediate to small x region relative to the nCTEQ15 results. While the shift of the PDF in the reweighting was in the same direction as in this analysis, its magnitude was much smaller.
We now turn our attention to the error band of the gluon distribution in Fig. 10. At NLO, the gluon enters for the first time the W ± and Z boson production through the gq initiated contributions. The addition of the W ± /Z LHC data to the fit is thus not expected to add significant constraining power for the gluon distribution. Contrary to this naive expectation, due to high center of mass energy and relatively small values of the probed x, the gluon distribution can have a considerable contribution to W ± /Z production processes; this is reflected in the reduced error bands of Norm3 as compared with nCTEQ15. Indeed, an independent variation of the open gluon parameters around the minimum in the Norm3 fit confirms that the χ 2 contribution from the LHC data is similarly steep or steeper than contribution from all the other data included in the fit. 9

A. Comparison with other nPDFs
Having investigated the impact of the W ± /Z heavy ion data including normalization effects, we now compare our PDFs with other results from the literature.
There are a number of nPDF sets available [60,65,66] including some new determinations [5,6,67]. The TUJU19 analysis [67] extends the xFitter framework to include nuclear PDFs; this open-source program provides a valuable tool for the PDF community. As an initial step, TUJU19 assumed s =s and s =ū =d, and the resulting nPDFs compare favorably with EPPS16 and nCTEQ15 within uncertainties.
A separate effort by the NNPDF collaboration [5,6] uses neural network techniques to extract the gluon and quark nPDFs; this method provides a complementary approach to the traditional parameterized function-based method. Their recent analysis [6] has produced the nNNPDF2.0 nPDF set which includes charged current DIS data from NuTeV (Fe) and Chorus (Pb), and also LHC W ± /Z data. They also compute the strangeness ratio, R s = (s +s)/(ū +d), and find the nuclear value is reduced as compared to the proton. The neutrino DIS data and LHC W + c associated production seem to prefer a lower R s value, while the inclusive W and Z production favor a larger value. These interesting observations raise some important issues, and additional investigation is warranted to better understand the strange distribution [7].
The EPPS16 data sets include DIS, DY, RHIC inclusive pion, and LHC W ± /Z and dijet data; in particular, this set incorporates a number of parameters to provide flexibility in both the strange and gluon PDFs. Therefore, it will be interesting to compare the variation of these flavors between our original nCTEQ15 nPDFs and our nCTEQ15WZ fit.
The nCTEQ15WZ fit is based on the Norm3 fit (with 3 normalization parameters), and in addition includes the RHIC pion data in the fitting loop. The RHIC pion data is fit with the Binnewies-Kniehl-Kramer (BKK) fragmentation functions [68] using a custom griding technique for fast evaluation [3]. The resulting nCTEQ15WZ nPDFs are nearly identical as the Norm3 nPDFs which is evident when comparing the χ 2 /N dof values of Table III, as well as the PDFs in Fig. 11.
We now compare the results of our nCTEQ15WZ fit with the nCTEQ15, EPPS16, and nNNPDF2.0 nPDFs in Figs. 12 and 13. To begin, we focus on the plots at Q = 2 GeV as the variations are more evident here. For the up and down components {u, d,ū,d}, nCTEQ15WZ is quite similar to nCTEQ15, and these flavors generally lie below EPPS16 and nNNPDF2.0, but are within uncertainties. For the strange and gluon, we see that nCTEQ15 and EPPS16 are generally similar for larger x values, and then diverge somewhat for small x. The nCTEQ15WZ nPDFs lie below nCTEQ15 and EPPS16 for large x values, and then above at intermediate to small x values; this allows s(x) and g(x) to increase the W ± /Z cross section in the region of the data (x ∼ 0.02) while not perturbing the momentum sum rules. nNNPDF2.0 is similar to nCTEQ15 and EPPS16 for large x values, but then increases for smaller x. For the strange distribution, nNNPDF2.0 coincides with nCTEQ15WZ at small x, while for the gluon, nNNPDF2.0 exceeds nCTEQ15WZ at small x. Similar effects to the above are generally evident at larger Q values (Fig. 13), but their magnitude is diminished due to the DGLAP evolution effects.

B. Comparison with proton results
The strange quark PDF has also been studied extensively for the proton case by many groups including ABM [11], CT18 [1], JAM [13], MMHT [14,15], and NNPDF [16]. There is a close connection between the proton and nuclear PDFs; for example, nCTEQ15 uses the proton PDF as a boundary condition, and EPPS16 fits nuclear ratios relative to the proton.
One quantity of interest we can compare between the proton and the nuclear PDFs is the ratio of the strange PDF relative to the light-sea quarks: R s = (s +s)/(ū +d). In Fig. 14, we compute R s for selected Q values, and compare this to the proton result as extracted by ATLAS [29,70].
Comparing the proton and the lead results at Q 2 = 1.9 GeV 2 , we see that the behavior of the Norm3 curve (panel-b) is quite similar to the proton result (panel-a). In contrast, the nCTEQ15 result is generally flat across all x values as the strange was set to be a fixed fraction of the u/d-sea PDFs, s =s = κ(ū +d)/2. Additionally, we also display the other fits, Norm0 and Norm2, to illustrate the range of possible variations. The uncertainty bands for Norm3 are displayed; these are large for small x, where the strange is poorly constrained, and also at very large x where the quark sea denominator vanishes. We also display larger Q values which illustrates the convergent effects of the DGLAP evolution.
In the previous section we raised the question as to whether the enhanced strange distribution was reflecting the true underlying physics, or was instead an artifact of the fit. The similarities of R s between the proton and lead PDFs may indicate that the enhanced strange PDF is, in fact, a real effect. To definitively answer this question will require additional analysis, and this work is ongoing.

V. CONCLUSION
Our ability to fully characterize fundamental observables, like the Higgs boson couplings and the W boson mass, and to constrain both SM and BSM signatures is strongly limited by how accurately we determine the underlying PDFs [71]. A precise determination of the strange PDF is an important step in advancing these measurements.
The new nCTEQ++ framework allowed us to include the LHC W/Z data directly in the fit. While these new fits significantly reduced the overall χ 2 for the W/Z LHC data, we still observe tensions in individual data sets which require further investigation. Our analysis has identified factors which might further reduce the apparent discrepancies including: increasing the strange PDF, modifying the nuclear correction, and adjusting the data normalization.
Compared to the nCTEQ15 PDFs, these new fits favor an increased strange and gluon distribution in the x region relevant for heavy ion W ± /Z production. While we obtain a good fit in terms of the overall χ 2 values, we must ask: i) how the uncertainties and data normalization affect the resulting PDFs, and ii) whether the results truly reflect the underlying physics, or is the fit simply exploiting s(x) because that is one of the least constrained flavors? The answer to this important question will require additional study; this is currently under investigation.

ACKNOWLEDGMENTS
We are pleased to thank Aaron Angerami,Émilien Chapon, Cynthia Keppel, Jorge Morfin, Pavel Nadolsky, Jeff Owens and Mark Sutton for help and useful  FIG. 14 (a) Recent (preliminary) result from ATLAS on the strange ratio for the proton [69]. (b) The nuclear strange ratio for lead (Pb) nPDFs as obtained in our fits. The uncertainty band for nCTEQ15 is shown in gray, and Norm3 in blue.

Appendix A: Fitting data normalizations
When fitting the normalization of data sets, we use the χ 2 prescription given in Ref. [72]. For a data set D with N data points and S correlated systematic errors, the χ 2 of the data set reads: where σ norm is the normalization uncertainty and T i is the theoretical prediction for point i. The last term of Eq.(A1) is called the normalization penalty and it enters when the fitted normalization, N norm , differs from unity. The normalization uncertainty σ norm appearing in the denominator prevents large excursions of N norm away from unity.
The covariance matrix C ij is defined as: where σ i is the total uncorrelated uncertainty (added in quadrature) for data point i, andσ iα is the correlated systematic uncertainty for data point i from source α. Using the analytical formula for the inverse of the correlation matrix as in Ref. [73], we obtain: and