Massive charged-current coefficient functions in deep-inelastic scattering at NNLO and impact on strange-quark distributions

We present details on calculation of next-to-next-to-leading order QCD corrections to massive charged-current coefficient functions in deep-inelastic scattering. Especially we focus on the application to charm-quark production in neutrino scattering on fixed target that can be measured via the dimuon final state. We construct a fast interface to the calculation so for any parton distributions the cross sections can be evaluated within milliseconds by using the pre-generated interpolation grids. We discuss agreements of various theoretical predicitons with the NuTeV and CCFR dimuon data and the impact of the results on determination of the strange-quark distributions.


Introduction
Deep-inelastic scattering (DIS) remains the most important probe of quantum chromodynamics (QCD) since it was pioneered almost fifty years ago. Various dedicated experiments have been carried out since then aiming to study the internal structure of nucleons, e.g., the HERA experiments by colliding electron or positron with proton. The HERA measurements on neutral-current (NC) and charged-current (CC) DIS provide the backbone constraint in modern determination of the parton distribution functions (PDFs) [1,2]. The heavy quark, especially charm quark plays an important role in describing the structure functions of proton measured in DIS within the framework of QCD factorization. In perturbative QCD calculations heavy-quark mass dependence of the DIS coefficient functions are crucial for analyses of the DIS data, in the inclusive structure function measurements and even more in the open production of heavy quarks, and ensure a precise determination of PDFs that is vital for the ongoing programs at the Large Hadron Collider (LHC). In the neutral-current case the DIS coefficient functions are known up to O(α 2 S ) with exact heavy-quark mass effects [3,4]. For the charged-current case the heavy-quark mass effects have been calculated to O(α S ) in Refs. [5][6][7], to approximate O(α 2 S ) in Ref. [8]. Recently the exact O(α 2 S ) results with full mass dependence have been completed [9]. The O(α 3 s ) results are also available for structure function xF 3 at large momentum transfer [10].
Among all DIS experiments there exist one specific measurement, charm-quark production in DIS of a neutrino from a heavy-nucleus. It provides direct access to the strange quark content of the nucleon which is poorly constrained by the inclusive structure function data. At lowest order, the relevant partonic process is neutrino interaction with a strange quark, νs → cX, mediated by the weak charged current. Experimentally one can require a semi-leptonic decay of the charm quark to muon that gives the so-called dimuon final state as measured by CCFR [11], NuTeV [12], CHORUS [13], and NOMAD [14] collaborations. In global determination of PDFs it is from those dimuon data which prefers a suppressed strangequark distribution than u and d sea-quarks [15,16,69]. That agrees with predictions from various models suggesting that the strange PDFs are suppressed compared to those of light sea quarks due to its larger mass [18][19][20]. The strange quark PDFs can play an important role in LHC phenomenology, contributing, for example, to the total PDF uncertainty in W or Z boson production [21,22], and to systematic uncertainties in precise measurements of the W boson mass and weak-mixing angle [23][24][25].
On another hand, thanks to the LHC we can also extract the strange quark PDFs independently from collider data only, e.g, via a combined analysis of HERA DIS data and the W , Z/γ * boson production data from the LHC. The latter can provide constraints on the strange quark PDFs due to its high precision and the fact that differential distributions can separate different sea flavors. The ATLAS collaboration have reported such a study, ATLAS-epWZ16 [26], using the HERA I and II combined data and the updated 7 TeV measurements on W and Z/γ * differential cross sections. Interestingly, an unsuppressed strange quark PDF is preferred with R s ≡ (s +s)/(ū +d) measured to be 1.13 +0.08 −0.13 at x = 0.023 and Q 2 = 1.9 GeV 2 . Similar conclusion has been reached in an earlier ATLAS study based on a smaller sample of W and Z data [27]. In comparison the values of R s from global PDF determination are 0.55 ± 0.21, 0.57 ± 0.17, 0.60 ± 0.13, and 0.63 ± 0.03 for CT14 [15], MMHT2014 [16], NNPDF3.1 [17], and ABMP16 [28]. In NNPDF3.1 they also performed an alternative fit using exactly the same data sets as in the ATLAS-epWZ16 analysis. They confirmed the pull on the central values of the strange-quark PDFs by the ATLAS data but arrived at a much larger uncertainties. The discrepancies seen between the determinations of strangeness from fixed-target dimuon data and the ATLAS data have attracted a lot of attentions recently. It was suggested in Ref. [29] that the ATLAS determination may be biased by the special parametrization form of PDFs adopted. While future LHC data might be helpful to clarify whether there are indeed tensions between those two determinations, it is also important to investigate various theoretical uncertainties especially in the case of charm quark production in DIS.
In Ref. [9] we have reported a first application of our O(α 2 S ) results on the massive coefficient functions of charged-current DIS. We calculated the next-to-next-to-leading-order (NNLO) QCD corrections to charm-quark production in DIS of a neutrino from a nucleon. The calculation is based on a phase-space slicing method and fully-differential Monte Carlo integration. We found the NNLO corrections can change the cross sections by up to 10% depending on the kinematic region considered. In this paper we provide further elaboration of the methods and numerical results of our NNLO calculation. Moreover, we implement our calculation into a fast interface based on grid interpolation that ensures repetition of the calculation within millisecond for arbitrary PDFs. It allows a first study of effects of the NNLO massive coefficient functions on determination of the strange-quark PDFs in the context of Hessian profiling, and can be used in future global analysis of PDFs.
In the remaining paragraphs we outline the method used in the calculation, present detailed numerical results on QCD corrections on cross sections of charm-quark production in charged-current DIS, demonstrate the accuracy of the grid interpolation, study the agreements of data and theory with various PDFs, and finally study effects of the NNLO corrections on extraction of the strange-quark PDFs.

NNLO calculation
We have presented briefly the framework of our NNLO calculation together with selected numerical results in Ref. [9]. Here we give more details on the theoretical ingredients of the calculation as well as more numerical results focusing on kinematic region of the fixed-target dimuon measurements. We also discuss the applicable kinematic range of our calculation utilizing fixed-flavor number scheme for heavy quarks and the possible improvement by extending to a variable-flavor number scheme.

Theoretical framework
The perturbative calculation utilizes a generalization of the phase-space slicing method to NNLO as motivated by the q T subtraction method proposed in [30]. There have been quite a few recent applications of similar methods on either decay processes [31], scattering at lepton colliders [32,33], and hadron colliders [34][35][36][37][38][39]. The key ingredient of above methods is the use of soft-collinear effective theory (SCET) and heavy-quark effective theory (HQET) [40][41][42][43] to systematically factorize the cross section and derive its perturbative expansion in fully unresolved region of QCD radiations as was proposed in Ref. [31]. Note here we adopt HQET for purpose of extracting the soft singularities in the perturbative expansion which is irrelevant to the actual mass of the charm quark. Other approaches on handling singularities in the fully unresolved region include sector-improved FKS subtraction method [44,45], sector decomposition method [46,47], antenna subtraction method [48,49], colorful subtraction method [50], and Projection-to-Born method [51]. Charm quark production at leading order (LO) through weak charged current in DIS can be represented by the diagram in Fig. 1. There are also Cabibbo suppressed contributions from d quark initial state. The production of charm anti-quark is similar. In our phase-space slicing method we first define a resolution variable which can isolate the unresolved phase space. As was discussed in Ref. [9], the appropriate resolution variable in this case is a fully inclusive version of beam thrust [52] or N-jettiness [53], which differs from the standard beam thrust or N-jettiness in that no partition in the phase space of final-state radiation is imposed, as there is only one collinear direction in the problem. In Eq. (2.1), p X is the momentum of total QCD radiation in the final state, p c is the momentum of the charm quark, and q is the momentum transfer as carried by a virtual W boson. p n is a momentum align with the incoming beam, whose large lightcone component equals the large lightcone component of the incoming momentum entering the W sc vertex. Here the lightcone direction n is chosen as the direction of the incoming beam, andn = (1, − n).
With the definition for τ , the differential cross section for any infrared-safe observable O can be separated into resolved and unresolved part, Further we can write down a factorization formula for the unresolved contribution, up to power corrections of the form τ cut ln k τ cut , where E d is the energy of the s quark entering the W cs vertex. The derivation of this factorization formula is very similar to the derivation of beam thrust of N-jettiness factorization [53]. In Eq. (2.3), dσ (0) (z)/dO is the Born level partonic differential cross section for the process 4) where P N is the momentum of the incoming hadron. The variable y is defined as y = q 2 /m 2 c < 0. The hard function H(y, µ) for charm quark production can be straightforwardly related to the hard function for bottom quark decay through analytic continuation. We refer readers to Refs. [54][55][56][57] for the full two-loop results 1 .
The soft function S is defined as a vacuum matrix element of Wilson loops. In a practical calculation, they can be obtained by taking the eikonal limit of the real corrections, with the insertion of a measurement function δ(k s − k·n), where k s is the total momentum of the soft radiation in the final state. For instance at one-loop the soft function can be calculated from the diagrams where the star distribution is defined as We refer to Ref. [58] for the full two-loop soft function. The beam function B is defined as the matrix element of collinear field in a hadron state (proton in our case), with the virtuality t = 2p n ·l of the beam jet measured [52], where l is the momentum of final state collinear radiation, and p n is defined in Eq. (2.1). The beam function can be written as convolution of perturbative coefficient functions and the usual PDFs, For example, the one-loop quark-to-quark coefficient function can be calculated through the diagrams (2.9) We also need the gluon-to-quark coefficient function at this order. The quark beam function has been calculated through to two loops [59]. We quote the result up to one-loop here (2.10) Substituting the expansion of hard, soft and beam function into the factorization formula in Eq. (2.3) gives the leading power in τ prediction for the unresolved distribution. The dependence on τ is very simple and can be integrated out analytically. Note that the power suppressed terms neglected in Eq. (2.3) also can be calculated analytically and used to improve convergence of the phase-space slicing method as demonstrated in Refs. [60,61]. For a small cut-off τ cut , integration of the unresolved distribution obtained from the factorization formula results in large logarithmic dependence on the cut-off. For sufficiently small cut-off, the large cut-off dependence is to be canceled by the resolved contribution, up to Monte-Carlo integration uncertainty. The resolved contribution, as its name suggests, is free of infrared singularities at NLO. At NNLO, the resolved contribution contains subdivergences. These sub-divergences cannot be resolved by our resolution variable τ . They must be canceled using other methods. Fortunately, the infrared structure of sub-divergences is lower by one order in α S than the unresolved part. For a NNLO calculation, we can use any existing subtraction method to cancel the sub-divergences. In our calculation, we employ the dipole subtraction formalism [62,63] to remove the sub-divergences. We also need the one-loop amplitudes for charm quark production with an additional parton, and tree-level amplitudes for charm quark production with two partons. We extract the former from Ref. [64]; for the later we use HELAS [65]. The calculations in resolved region have been cross checked with Gosam [66] and Sherpa [67] and full agreement are found.

Numerical results
We move to numerical results for the reduced cross sections of charm-quark production in DIS of neutrino on iron. We use CT14 NNLO PDFs [15] with n f = 3 active quark flavors and the associated strong coupling constant by default. We use a pole mass m c = 1.4 GeV for the charm quark, and CKM matrix elements |V cs | = 0.975 and |V cd | = 0.222 [68]. The renormalization and factorization scales are set to µ 0 = Q 2 + m 2 c unless otherwise specified. We choose a phase-space slicing parameter of τ cut = 10 −3 which is found to be small enough to neglect the power corrections [9].
In Fig. 2 we show the QCD corrections to a differential reduced cross section in Bjorken x for which the electroweak couplings have been taken out. We plot the NLO and NNLO The scale variations at LO are large in general but vanish at x ∼ 0.1 indicating its limitation as estimate of perturbative uncertainties. It was found even the scale variations at NLO underestimate the perturbative uncertainties at small and moderate x regions due to accidental cancellations as will be explained later. The NNLO scale variations give a more reliable estimation of the perturbative uncertainties and also show improvement at high-x compared with the NLO case. Results are similar for charm anti-quark production which can be related via a charge conjugate parity transformation except for the differences of initial state PDFs. Especially the charm quark production involves Cabibbo suppressed contributions at tree level from d-valence quark which dominate at high-x while only sea-quark contributions exist for charm anti-quark production.
In the small and moderate x region the NNLO corrections are almost as large as the NLO corrections. That motivates a careful examination of the convergence of the perturbative expansion. In Fig. 3 we plot the QCD corrections from two main partonic channels, i.e., with the strange (anti-)quark initial state, including Cabibbo suppressed d(d) quark contributions, and with the gluon initial state, for the same distribution as shown in Fig. 2. The right plot of Fig. 3 shows the corrections for charm anti-quark production. We observe a strong cancellation among the NLO corrections from the strange anti-quark and the gluon channels starting from small-x and persisting to high-x region. We regard this cancellation accidental in that it does not arise from basic principles but is a result of several factors. The cancellation of NLO corrections remains if instead using the NLO PDFs or alternative NNLO PDFs e.g., Figure 3: QCD corrections at different orders separated into partonic channels for a differential reduced cross section in Bjorken x for charm (anti-)quark production from (anti-)neutrino scattering on iron target.
MMHT2014 [16] and NNPDF3.0 [69]. A similar cancellation has also been observed in the calculation for t-channel single top quark production [70]. The size of NNLO corrections are smaller than the NLO ones for the individual partonic channels indicating good convergence of the perturbative expansion. However, the cancellation between the two channels is much mild at NNLO which results in a net correction as large as the NLO one. For this reason we expect the corrections from even higher orders to be smaller than the NNLO corrections. The left plot in Fig. 3 shows results for charm quark production for which the situation is similar at low-x region. At high-x the correction from gluon channel flattens out due to the smaller size of gluon PDF as comparing to d-valence PDF, and the net correction is small and positive at NNLO. We further calculate the double differential reduced cross sections in x and y as was measured by various experimental groups. We choose the kinematics and neutrino energies as those in CCFR [11] measurement. In Fig. 4 we plot ratios of various predictions to the LO differential cross sections for charm quark production with three different energies and each with three choices of y. Here we use a charm-quark mass of 1.3 GeV. The solid and dotted curves correspond to using scales of µ 0 and 2µ 0 . Note the LO cross sections in the denominator are always evaluated with the scale µ 0 . For the nominal scale choice (µ 0 ) the NNLO corrections are about −10% at x ∼ 0.02 and a couple of percents at x ∼ 0.3. The size of QCD corrections increases with y in low-x regions. Dependence on the beam energy is in the opposite direction and is weaker in general. Scale dependence of the NNLO predictions are slightly weaker than those of the NLO predictions at small-x. In moderate and large-x regions the NLO predictions show a scale dependence that is too small due to the strong cancellations mentioned earlier. Fig. 5 shows similar results for charm anti-quark production. The QCD corrections are even more pronounced in this case due to the relatively larger gluon contributions. For y = 0.802 the NNLO corrections can reach −15% for x ∼ 0.02 and remain −10% for x ∼ 0.2. The same conclusion holds for the scale dependence as in the case of charm quark production. Fe c+X with CT14NNLO_NF3 Figure 4: QCD predictions at different orders with scale choices of µ 0 and 2µ 0 for a double differential reduced cross section in Bjorken x and inelasticity y for charm quark production from neutrino scattering on iron target.
Since the charm quark production is usually measured at low to moderate momentum transfers, the theoretical predictions can depend significantly on the choice of charm-quark mass, in our case the pole mass. Note the determination of charm-quark pole mass has an intrinsic uncertainty of 0.1 ∼ 0.2 GeV due to the renormalon ambiguity. In Fig. 6 we show the ratio of double differential cross sections calculated when using a charm-quark pole mass of 1.5 GeV to 1.3 GeV, at LO, NLO, and NNLO, for charm anti-quark production. The results for charm quark production are similar and not shown for simplicity. At LO the charmquark mass dependence can be calculated easily. The dominant part of that is known as slow rescaling [71] due to the kinematic suppression, i.e., by replacing the momentum fraction in evaluation of PDFs with ξ = x(1 + m 2 c /Q 2 ). That explains the trends we show in Fig. 6. The cross sections with larger charm-quark mass are especially suppressed in small-x region and for smaller neutrino energies where the Q 2 is low. Shapes of the suppression factor with respect to x are different for different values of y due to the sub-dominant dependence on the mass from the hard matrix elements. The mass dependence is insensitive to higher order corrections. The NLO predictions show a slightly weaker suppression comparing to LO ones Fe c+X with CT14NNLO_NF3 Figure 5: Similar as Fig. 4 for charm anti-quark production from anti-neutrino scattering.
in general, especially at large-x and smaller neutrino energies. Effects of NNLO corrections on the mass dependence are almost negligible for the full kinematic range considered.

Heavy quark scheme
The NNLO calculations are carried out in a fixed flavor number scheme with n f = 3. This should be the appropriate scheme for Q m c . For the semiinclusive charged-current (CC) DIS process we studied, at Q m c , there exist logarithmic contributions of ∼ α n S ln n (Q 2 /m 2 c ) due to initial-state gluon splitting into a cc pair in the quasi-collinear limit. 2 In principle that needs to be resummed by using the heavy-quark PDFs together with an appropriate general mass variable flavor number (GM-VFN) scheme, for example, the ACOT [72], FONLL [73], or RT [74] schemes. The exact O(α 2 S ) massive coefficient functions [9] complete all ingredients needed for constructing such a scheme like S-ACOT-χ [75] at NNLO for the charged-current scattering. For the kinematics where the dimuon measurements carried out, the Q 2 is not too high comparing to the charm-quark mass in the bulk of the data. Besides, the experimental uncertainties are at least at the level of 5 ∼ 10% for the NuTeV and CCFR measurements. Fe c+X with CT14NNLO_NF3 Figure 6: Dependence of a double differential reduced cross section in Bjorken x and inelasticity y on the charm quark mass, shown in ratios of predictions with m c = 1.5 GeV to 1.3 GeV, for charm quark production from neutrino scattering on iron target.
Thus a VFN scheme is not of immediate relevance for the phenomenological study on dimuon measurements. We leave a formal study on the GM-VFN scheme for future publications while providing an estimate of the logarithmic contributions beyond O(α 2 S ) in below. As mentioned earlier the logarithmic contributions can be resummed effectively with a perturbative charm (anti-)quark PDF in n f = 4 scheme that follows a DGLAP evolution, where P ij is the DGLAP splitting function with dependence on n f suppressed, µ is the factorization scale. The exact results for P ij are known up to three loops [76,77]. The charmquark PDF at arbitrary scales can be derived from the boundary conditions at µ = m c by evolving upward. Note that starting at O(α 2 S ) the charm-quark PDF has a small discontinuity at µ = m c . We can expand the charm-quark PDF in the strong coupling constant, where ∆ (2) is of O(α 2 S ) due to the discontinuity crossing the heavy-quark threshold and L = ln(µ 2 /m 2 c ). It is understood that the strong coupling constant, the one-and two-loop splitting functions P (0,1) ij , and the one-loop β function in Eq. (2.12) are all evaluated with n f = 4. We can translate the strong coupling constant and PDFs with n f = 4 to those with n f = 3 via matching at the charm-quark threshold [78]. Furthermore, we can expand them in α S (µ 2 ) instead. We arrive at an expanded solution, with β functions and splitting functions for n f = 4. Those O(α S ) and O(α 2 S ) logarithmic contributions have already been captured by our NLO and NNLO calculations respectively. Differences of the evolved charm-quark PDF and the expansion in Eq. (2.13) can serve as an estimate of the remaining logarithmic contributions at higher orders.
In Fig. 7 we plot the differences of an evolved charm-quark PDF c evol. at NNLO (with 3loop splitting functions) and the expanded solution c exp. up to O(α S ) and O(α 2 S ) as a function of µ = Q for several x values, normalized to the effective strangeness PDF s ef f which is a combination of s(s) and d(d) PDFs. The charm quark production cross section at LO is simply proportional to the effective strangeness PDF with the slow rescaling. We use the MSTW2008 NNLO PDFs [79] with m c = 1.4 GeV as an input. For small x values we can see the FFN calculation at O(α S ) misses a large portion of the logarithmic contributions that can reach 10-20% of the LO charm quark production cross sections for Q ∼ 10 GeV. On another hand the NNLO calculation can reproduce well the resummed contributions with the remaining logarithmic contributions of about 2% of the LO cross sections for the same Q values. Note the highest Q value that the CCFR and NuTeV measurements probed is around 10 GeV. For large x values the conclusion is similar for charm quark production. For charm anti-quark production the charm-quark PDF has a relatively larger weight due to the quick falling of the sea-quark PDFs. The contributions beyond NNLO can reach 5% for x = 0.3 and Q = 10 GeV.

Fast interface
The above calculation can not be immediately used in the global analysis of QCD due to the time-consuming nature of NNLO calculations and the fact that the analysis involves scan over a large number of PDF ensembles. Indeed even the NLO calculation is inadequate for direct use in the analysis. The PDF fitting group instead needs to rely on either K-factor approximation or fast interface based on grid interpolations. There have been quite some developments on fast interface for high-order perturbative calculations, e.g., APPLgrid [80], FastNLO [81], aMCfast [82], starting from NLO in QCD to NNLO most recently [83]. We have constructed a fast interface specialized for our calculation following similar approaches. First of all the PDFs at arbitrary scales can be approximated by an interpolation on a onedimensional grid of x, where we choose the interpolation variable y(x) = x 0.3 and the interpolation order n = 4, and f j is the PDF value on the j-th grid point. δy has been chosen so as to give 50 grid points between x = 1 and a minimum determined according to the specified kinematics. We use a n-th order polynomial interpolating function I (n) i and the starting grid point k is determined so as x located in between the (k + 1)-th and (k + 2)-th grid points. The cross section in deep inelastic scattering can thus be expressed as where the summation runs over different sub-channels p, perturbative orders m, and the grid points i. The interpolation coefficients B(p, m, i) which are independent of the PDFs can be obtained by projecting the event weight onto the corresponding grid points during the MC integration. Once those interpolation coefficients are calculated and stored, the cross sections with any PDFs can be obtained via Eq. (3.2) without repeating the time-consuming calculations of the matrix elements.
In Table 1 we show the typical time cost in a direct calculation and the interpolation with the NuTeV kinematics. Also shown is the time cost for generating the interpolation grid. The direct calculation involves intensive MC integration as expected, costing about 60 CPU core-hours per data point of the double differential distribution d 2 σ/dxdy in charmquark production with NuTeV kinematics. The grid generation costs four times more since it requires separation of different sub-channels. However, with the generated grid, for any PDFs the interpolation/calculation takes less than a millisecond. The precision of the interpolation are found to be around a few permille at NNLO, smaller than the typical errors from MC integration. In Fig. 8 the solid line shows the ratio of the cross sections from direct calculation and the fast interpolation using the grid generated from the same run both using CT14 NNLO PDFs, for all the data points in NuTeV and CCFR measurements with charm (anti-)quark production. In this case deviation of the two predictions are simply due to the interpolation errors. Also shown in Fig. 8 are comparison of the interpolation results for MMHT2014 and NNPDF3.0 NNLO PDFs with the grid generated from CT14, with independent direct calculations using the same PDFs. Here the two predictions for each PDF choice differ at most half percent due to the MC integration errors in the direct calculations as shown by the error bars.

Impact on strange-quark distributions
Now we move to discuss potential impact of the NNLO calculations on constraining parton distributions, especially the strange-quark distributions, by checking the agreements of different theories with experimental data. We select the NuTeV and CCFR measurements on charm (anti-)quark production in the form of double differential cross sections d 2 σ/dxdy, which provides dominant constraints on the strange-quark distributions, e.g., in MMHT2014 and CT14 global analyses. The theoretical predictions used in those analyses are at NLO only. We include the data point only with Q 2 > 4 GeV 2 to leave out the region where higher-twist corrections can be potentially large. That results in 38(33) data points for charm (anti-)quark production in NuTeV, and 40(38) data points for charm (anti-)quark production in CCFR. Besides, we have simply corrected the data for nuclear effects of iron target [22,84] using a parametrization of F 2 ratio measured at SLAC and NMC, instead of including more sophisticated corrections to individual parton flavors [85][86][87][88] in the theory calculations. That leads to corrections on the data of 2% at x ∼ 0.05, -4% at x ∼ 0.1, and 5% at x ∼ 0.4. We did not include uncertainties on the nuclear corrections since the correction itself is already small comparing to experimental errors for the x range considered. The experimental uncertainties include the total statistical and systematic errors which are treated uncorrelated among different data point. The total error for each data point has been scaled by square root of its effective freedom so as a reasonable fit should have χ 2 /N pt of one [12]. Besides, there is an additional systematic error of 10% assumed due to the input of semi-leptonic decay branching ratio of charm quark when unfolding the dimuon cross sections back to the charm production cross sections [12]. This normalization error is assumed to be fully correlated among all data points in both NuTeV and CCFR measurements. In the following we first compare predictions with various PDFs to the experimental measurement. Later we show how the strange quark distributions may change if using our NNLO results instead of the NLO ones in the PDF analyses by means of Hessian profiling [89].

Theory-data agreement with various PDFs
We considered the updated NNLO PDFs from the major PDF fitting groups, including CT14 [15], MMHT14 [16], NNPDF3.1 (both nominal set and set with only collider data) [17], ABMP16 [90], HERAPDF2.0 [1], and ATLAS-epWZ16 [26]. For the case of NNPDF, the original PDF representation of MC replicas has been transformed into a Hessian PDF set using the MC2H package [91]. Note that we have used NNLO PDFs for both NLO and NNLO calculations. In case that the PDFs with n f = 3 are not publicly available we evolve the nominal PDFs by DGLAP evolution at 3-loop and with n f = 3 starting from a scale below the charm quark mass threshold.
We show the fits to NuTeV and CCFR data (149 data points in total) with NLO and NNLO predictions for various choices of PDFs in Table 2. Here through the calculations we have used a charm-quark pole mass of 1.3 GeV and a scale of µ = Q 2 + m 2 c despite the fact that different PDF groups use a charm-quark pole mass ranging from 1.3 to 1.6 GeV. 3 In each fit we show the χ 2 and the normalization shift (in unit of 1 σ error) of the central PDFs without and with including the full PDF uncertainties. In the later case it gauges the overall agreement between data and theory with both uncertainties included. Shift with minus sign indicates the data prefer smaller values for the cross sections. Each pair of eigenvector PDFs corresponds to one correlated theoretical error with symmetric Gaussian distribution when including the PDF uncertainties [89]. For the HERA and ATLAS fits the PDF uncertainties include those model and parametrization uncertainties as well. Note the χ 2 shown here may not represent the actual fit quality of the same data in their global analyses since there different predictions or input parameters are used.
The PDFs shown in Table 2 fall into two distinct groups, those without including any dimuon data in the PDF analysis, namely the HERA, ATLAS, and NNPDF collider-only PDFs, and the others with dimuon data, either from NuTeV, CCFR, CHORUS, or NOMAD. The χ 2 /N pt are about one for PDFs in the second group as all the dimuon data consistently prefer a suppressed strangeness as discussed in the introduction section. Interestingly, CT14, MMHT14 and NNPDF3.1 show very similar results on χ 2 and also the normalization shift. The NNLO predictions without PDF uncertainties give slightly smaller χ 2 in general for the specified mass and scale as comparing to NLO, though those PDFs are fitted with NLO or approximate NNLO predictions. For PDFs from collider-only data, the fits are rather poor if  Table 2: χ 2 and normalization shift (in unit of 1 σ error) of fits to NuTeV and CCFR charm production data with various theoretical predictions using m c = 1.3 GeV and µ = Q 2 + m 2 c . The shifts are shown in brackets with minus sign indicating the data prefer smaller values for the cross sections. The numbers in bold font correspond to fits including the full PDF uncertainties as well.
not taking into account the PDF uncertainties. One reason is the ATLAS W/Z data do prefer larger central values of the strange-quark PDFs. With the PDF uncertainties included, the χ 2 /N pt have been reduced to below one for HERA and NNPDF collider-only PDFs indicating consistency of their PDFs and the dimuon data once both uncertainties are considered. The situation is different for the ATLAS PDFs where χ 2 /N pt is still about 1.5 even for the NNLO predictions. That can be further visualized by a direct comparison of theory and data as in Figs. 9 and 10 for NuTeV charm quark and anti-quark production respectively. Most of the data points lie far outside the PDF error bands with a non-trivial shape dependence. The PDF uncertainties of charm quark cross sections are smaller than the ones of charm anti-quark in general since the former also involves contribution from d valence quark which is better constrained than sea quarks. We conclude the ATLAS PDFs can not describe the dimuon data well and the NNLO calculations can only bring in limited improvement.
We further compare predictions from the same PDF sets but with different scale and charm quark mass inputs as shown in Tables 3 and 4. From Table 3 we can see that by using a scale of twice of the nominal choice the agreement between NLO predictions and data deteriorates, especially in the case of without PDF uncertainties. In comparison the χ 2 for NNLO predictions are less sensitive to the change of scale for both with and without PDF uncertainties. The cross sections are reduced especially at small-x when using a larger charmquark mass of 1.4 GeV as shown already in Fig. 6. As show in Table 4 that leads to a smaller χ 2 at NLO in general comparing with Table 2 when no PDF uncertainties are included. At NNLO the χ 2 can either decrease or increase depending on the PDFs considered indicating a different preference of charm-quark mass at NNLO comparing with NLO for certain PDFs. The χ 2 for the ATLAS PDFs are reduces as well. However, in both cases the χ 2 are still over 200 for predictions with the ATLAS PDFs.    Table 4: Similar as Table 2 but with m c = 1.4 GeV and µ = Q.

Hessian profiling of PDFs
One main motivation of this paper is to investigate impact of the NNLO calculations on extraction of the strange quark PDFs in global analyses including the dimuon data. Precisely we would like to see how the outcome strange quark PDFs change when using NNLO predictions instead of the NLO ones, as only NLO predictions are available in previous PDF fits. That could be done by individual PDF groups using the fast interface and grids presented in this paper. Alternatively we can estimate the possible shift of the PDFs by means of Hessian profiling [89]. In Hessian profiling the PDF parameters are assumed to follow a prior multi-Gaussian distribution. That corresponds effectively to a parabolic shape of the prior χ 2 around the central PDF and with ∆χ 2 = 1 when reaching the 1 σ error in each eigenvector direction. The χ 2 of any new data set to be included also forms a parabola in the PDF parameter space under linear approximation. The profiled PDFs are thus determined according to profiles of the total χ 2 , e.g., by minimization of the total χ 2 for the central value and ∆χ 2 = 1 for the PDF uncertainties. We stress here the Hessian profiling can only serve as an estimation of the effects of new data or theory on the PDFs since an actual PDF fit involves further complexities due to e.g., parametrization dependence and the requirement of tolerance conditions. We start with the HERAPDF2.0 NNLO PDFs which does not include any dimuon data but rather implements certain model constraints on the strangeness fraction and shape. We show the PDF ratio of strange to non-strange sea-quarks, R s at a scale of 2 GeV in Fig. 11. We can see a moderate suppression of the strangeness in small to intermediate x region comparing to the u and d sea-quarks and a rapid falloff at large x. The PDF uncertainties as indicated by the solid black lines are large, more than 30% in the entire x range, which include the experimental, model, and parametrization uncertainties. The colored bands in Fig. 11 are for the profiled PDFs with the NuTeV and CCFR charm (anti-)quark production data together using the NLO and NNLO theoretical predictions. It is clear that the dimuon data prefers an even suppressed strangeness with R s of about 0.5 in the full range of x. The profiled PDFs lie at the lower edge of the 1 σ error of the original PDFs indicating reasonable agreement between original PDFs and the dimuon data as already seen in Table 2. The profiled PDFs have a much smaller uncertainties on R s than the original PDFs as one expect. We notice that the PDF uncertainties are also reduced significantly in the small x region 10 −4 − 10 −2 which are beyond the coverage of the dimuon data. That is possibly due to the restricted parametrization form of strange quark PDFs used in the HERA PDF analysis. Importantly we found the NNLO predictions prefer higher values of R s than the NLO ones, in this case well above the 1 σ error band of the later. That can be understood since the NNLO corrections are negative in the bulk of the data. MMHT2014nnlo68cl_nf3 (1 PDF unc. ) Figure 12: Similar as Fig. 11 for profiling of the MMHT2014 NNLO PDFs.
We perform another profiling study with the MMHT2014 NNLO PDFs as shown in Fig. 12. Noted that since the MMHT2014 analysis already includes above dimuon data, the study here only means for checking impact of the NNLO corrections. We can see the NNLO predictions prefer larger strangeness than NLO predictions for x up to a few times 0.1 and by a similar amount as in Fig. 11. The shift of central values of the NLO profiled PDFs comparing to the original PDFs, though still within the PDF error band, is due to several facts. In the MMHT2014 fits [16] they use a charm-quark pole mass of 1.4 GeV and a semi-leptonic decay branching ratio of charm quark that is 7% lower than the one extracted by NuTeV and CCFR, both of which lead to an increase of the strange-quark PDFs. Besides, there are also LHC data in the MMHT analysis that pull the ratio further up. The uncertainties are largely reduced in the profiled PDFs mostly because we use the ∆χ 2 = 1 criterion rather than a dynamic tolerance condition as in the MMHT analysis. We have also compared the profiled PDFs with alternative scale choices and found those with NNLO predictions are less sensitive to the choice of scale.

Conclusion
In conclusion we have presented details on calculation of next-to-next-to-leading order QCD corrections to massive charged-current coefficient functions in deep-inelastic scattering. We focus on the application to charm-quark production in neutrino scattering on fixed target that can be measured via the dimuon final state as in the NuTeV and CCFR experiments. We construct a fast interface to the calculation so for any parton distributions the dimuon cross sections can be evaluated within milliseconds by using the pre-generated interpolation grids. The NNLO predictions thus can be conveniently included in future global analyses of QCD involving the dimuon data. We further compare the dimuon data with the NNLO predictions using various PDFs and confirm the pull of the ATLAS data on the strange-quark PDFs. Moreover, we study impact of the NNLO corrections on the extraction of strangequark PDFs in the context of Hessian profiling, and find with the NNLO predictions the dimuon data tend to favor larger strange-quark PDFs than with the NLO predictions. The fast interface together with the interpolation grids for dimuon cross sections are publicly available upon request. A definite conclusion on the potential inconsistency of the DIS and ATLAS data awaits the ongoing global fits by the PDF analysis groups.