Precision determination of the strong coupling constant within a global PDF analysis

We present a determination of the strong coupling constant \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha _s(m_Z)$$\end{document}αs(mZ) based on the NNPDF3.1 determination of parton distributions, which for the first time includes constraints from jet production, top-quark pair differential distributions, and the Z \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_T$$\end{document}pT distributions using exact NNLO theory. Our result is based on a novel extension of the NNPDF methodology – the correlated replica method – which allows for a simultaneous determination of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha _s$$\end{document}αs and the PDFs with all correlations between them fully taken into account. We study in detail all relevant sources of experimental, methodological and theoretical uncertainty. At NNLO we find \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha _s(m_Z) = 0.1185 \pm 0.0005^\text {(exp)}\pm 0.0001^\text {(meth)}$$\end{document}αs(mZ)=0.1185±0.0005(exp)±0.0001(meth), showing that methodological uncertainties are negligible. We conservatively estimate the theoretical uncertainty due to missing higher order QCD corrections (N\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^3$$\end{document}3LO and beyond) from half the shift between the NLO and NNLO \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha _s$$\end{document}αs values, finding \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta \alpha ^\mathrm{th}_s =0.0011$$\end{document}Δαsth=0.0011.


Introduction
The value of the strong coupling constant α s (m Z ) is a dominant source of uncertainty in the computation of several LHC processes. This uncertainty is often combined with that on parton distributions (PDFs), with which it is strongly correlated. However, while PDF uncertainties have reduced considerably over the years, as it is clear for example by comparing the 2012 [1] and 2015 [2] PDF4LHC recommendations, the uncertainty on the α s PDG average [3] remains substantially unchanged since 2010 [4]. As a consequence, the uncertainty on α s is now the dominant source of uncertainty for several Higgs boson production cross-sections [5].
Possibly the cleanest [6,7] determinations of α s come from processes that do not require a knowledge of the PDFs, such as the global electroweak fit [8]. These are free from the need to control all sources of bias which may affect the PDF determination and contaminate the resulting α s value. A determination of α s jointly with the PDFs, however, has the advantage that it is driven by the combination of a large number of experimental measurements from several different processes. This is advantageous because possible sources of uncertainties related to specific measurements, either of theoretical or experimental origin, are mostly uncorrelated amongst each other and will average out to some extent in the final α s result. In addition to the above, the simultaneous global fit of α s and the PDFs is likely to be more precise and possibly also more accurate than individual determinations based on pre-existing PDF sets, many of which have recently appeared [9][10][11][12][13][14][15]. This is due to the fact that it fully exploits the information contained in the global dataset while accounting for the correlation of α s with the underlying PDFs.
Here we present a determination of α s which exploits the most recent PDFs obtained with the NNPDF methodology, namely NNPDF3.1 [16]. This updates a previous determination of α s [17,18] based on NNPDF2.1 [19,20]. In comparison to this previous PDF set, NNPDF3.1 represents a substantial improvement both in terms of input dataset, theoretical calculations, and fitting methodology. Specifically, NNPDF3.1 is the first PDF set to make such an extensive use of LHC data as to be dominated by them. It is in fact the first global analysis to simultaneously use differential top, inclusive jet, and Z p T distribution data, all using exact NNLO theory. Indeed, typical PDF uncertainties are of order of one to three percent in the data region for NNPDF3.1, about a factor two smaller than they were for NNPDF2. 1. This greater precision in the PDF determination requires a corresponding improvement in the methodology used for the α s extraction. In our previous work [17,18], PDF replicas were determined for a number of fixed values of α s , which was then extracted from the χ 2 profile versus α s of the best fit PDF, obtained as an average over the replicas. Here instead, both α s and PDFs are determined from a simultaneous minimization in their combined parameter space. As we will discuss below, this new method corresponds roughly to determining the value and uncertainty on α s from the error ellipse of the multivariate measurement in the (α s , PDF) hyperspace, and the old method corresponds to performing a scan of α s along the best-fit PDF line, see Fig. 1 for a Fig. 1 Comparison between the standard deviation of a pair of correlated variables (α s , θ) and the one-sigma range for the variable α s along the best-fit line of θ. The best fit is denoted as (α s ,θ) and the ellipse is the one-sigma contour about it. The standard deviations on (α s , θ) are (σ α , σ θ ), while σ old is the one-sigma interval for α s with fixed θ =θ schematic illustration. In a situation when the variables are highly correlated, especially if the semi-axes of the ellipse are of very different length, the procedure used in our previous work might lead to an underestimate of the uncertainty in α s . Hence the new procedure becomes very relevant, now that some PDF uncertainties are rather small.
It turns out that the implementation of this simultaneous minimization within the NNPDF methodology is nontrivial: it requires the development of a suitable generalization of the standard NNPDF approach, which we call the correlated replica method. Using this strategy, α s can be treated like any other quantity that depends on the PDFs. In particular, its central value and uncertainty can be determined by performing statistics over a replica sample. This means that, for example, the uncertainty on α s is the standard deviation of an ensemble of α s values. As we shall see, this allows for a determination of α s with small experimental uncertainties, and negligible methodological uncertainties. Having reduced very much the size of all other uncertainties, the problem of accurately estimating theoretical uncertainties becomes quite serious. This is specifically problematic in the case of missing higher-order uncertainties (MHOUs), for which no fully satisfactory method has been developed. Here we will conservatively estimate the theoretical uncertainty due to missing higher order QCD corrections (N 3 LO and beyond) from half the shift between the NLO and NNLO α s values.
This paper consists of two main parts. First, in Sect. 2 we present the correlated replica method used for the determination of α s , explain how it is used to estimate the associated PDF uncertainties, and compare it with the method used in previous NNPDF determinations. Then, in Sect. 3 we present our determination of α s at NLO and NNLO together with a careful assessment of all sources of uncertainty. Possible future developments are briefly outlined in Sect. 4.

The correlated Monte Carlo replica method
As discussed in the introduction, the α s determination presented here differs from our previous one [17,18] because now the value of α s and its uncertainty are determined from a correlated fit together with the PDFs. After briefly summarizing the main aspects of the NNPDF methodology and the way it was used to determine α s in Ref. [17,18], we describe the main idea of the new method, and then discuss the details of its implementation. Only the salient aspects of the NNPDF methodology will be recalled here; the reader is referred to the original literature (see Ref. [16], of which we follow the notation, and references therein) and recent reviews [2,21,22] for a more detailed discussion.

General strategy
The NNPDF fitting methodology is based on constructing a Monte Carlo representation of the original data sample consisting of pseudodata (Monte Carlo replicas of the original data), and fitting PDF replicas to these data replicas. Specifically, starting with an N dat -component vector of experimental points D with components D i , a set of N rep replicas D (k) of the data is generated by means of: and σ stat i are normalization, systematic and statistical uncertainties, and r i are random numbers such that statistics over the replica sample reproduces the original statistical properties of the data in the limit of large N rep . For example, this means that where cov denotes the covariance over the replica sample and C i j is the full experimental covariance matrix of the data. A PDF replica is then fitted to each data replica D (k) . In the NNPDF approach, PDFs are parametrized using neural networks, in turn specified by a vector of parameters θ . In the most recent NNPDF3.1 analysis, this vector θ has 296 components, corresponding to 37 parameters for eight neural networks (for the up, antiup, down, antidown, strange, antistrange, total charm and gluon PDFs). Thus, for each data replica D (k) a best-fit θ (k) is found by minimizing a figure of merit characterizing the agreement between theory and data: Here, T i [θ ] is the theoretical prediction for the ith datapoint, dependent on the set of parameters θ , and C t 0 is the covariance matrix used in the fit. Recall that in the presence of multiplicative uncertainties, C t 0 cannot be directly identified with the experimental covariance matrix C used for pseudodata generation Eq. (2.1) lest the fit be biased [23], and must thus be constructed instead using a suitable procedure such as the t 0 method [24] (see also [25]). A peculiarity of the NNPDF approach is that the best-fit parameters of each replica, θ (k) , are not defined as the absolute minimum of the χ 2 Eq. (2.3) in order to avoid overfitting, i.e. in order not to fit statistical fluctuations. Instead, a suitable cross-validation algorithm is employed [26]. We thus obtain a set of best-fit PDF replicas θ (k) , each determined as the minimum with respect to θ of the figure of merit χ 2(k) computed using the kth data replica: (2.4) where argmin should be understood as minimization through cross-validation, rather than as the absolute minimum. Note that, because we employ non-deterministic minimization algorithms, specifically genetic algorithms, the best-fit θ (k) corresponding to a given data replica D (k) is not unique; two identical data replicas D (k 1 ) = D (k 2 ) may lead to two different θ (k 1 ) = θ (k 2 ) in two runs of the minimization algorithm.
In summary, the standard NNPDF methodology produces a set of replicas D (k) of the original data, and uses them to construct a set of PDF replicas which correspond to parameters θ (k) , where k runs over the replica sample.
The theory predictions T i , which enter in the figure of merit of the fit Eq. (2.3) depend not only on the PDF parameters θ , but also on theory parameters, specifically the value of α s . Therefore, in general we can view the figure of merit as a function χ 2 (α s , θ, D). In standard NNPDF determinations, α s is treated as a fixed parameter, along with all other theory parameters, such as quark masses, CKM matrix elements, the fine structure constant, and so on. On the other hand, it is well known (see e.g. Ref. [27] for an early reference) that PDFs are strongly correlated to the value of α s , so a determination of the combined PDF+α s uncertainty on a process which depends on both, requires knowledge of the PDFs as α s is varied. With this motivation, NNPDF sets are routinely released for different fixed values of α s , where the procedure of generating data replicas D (k) and determining PDF replicas determined by the best-fit parameters θ (k) is repeated several times for different values of α s .
In our previous work [17,18], α s was determined by first producing PDF fits for a range of values of α s . The χ 2 (α s ) of the mean of all the replicas was then fitted to a parabola as a function of α s . This methodology has two main drawbacks. The first is that, as mentioned, the PDFs are strongly correlated to the value of α s . With this method, however, the χ 2 profile is determined as a function of α s along the line in θ space which corresponds to the best-fit θ at each particular value of α s , without taking into account the variations in θ space. Hence, as illustrated in Fig. 1, with the methodology of Refs. [17,18] the resulting uncertainty on α s could be somewhat underestimated.
The second drawback is more subtle. In the NNPDF procedure, the PDF uncertainty is determined from statistics over the replica sample, so a one-sigma interval is determined by computing a standard deviation over replicas. Whether or not this corresponds exactly to a one-sigma (i.e. χ 2 = 1) interval in α s space is unclear. In fact, in PDF determinations based on Hessian minimization in parameter space, the χ 2 = 1 criterion is modified by a suitable tolerance factor [28][29][30], which possibly accounts for data inconsistencies or parametrization bias. It is unclear, but certainly possible, that PDF uncertainties estimated in the NNPDF fits also include, at least to some extent, such a tolerance.
Ideally, we would like a method of determining α s in which the uncertainty on α s is determined on exactly the same footing as the PDF uncertainty, and which thus yields the full probability distribution for α s , marginalised with respect to the PDF parameters. The goal is to treat α s on the same footing as the vector of parameters θ that determine the PDFs, i.e. to simultaneously minimize the figure of merit with respect to both α s and θ . This is difficult in practice, because the dependence on α s appears in the theoretical predictions, which, for reasons of computational efficiency, are provided in the form of pre-computed grids determined before the fit using the APFELgrid framework [31,32].
This difficulty can be overcome through the correlated replica method, as we now explain. The method relies on the concept of "correlated replica", or c-replica for short. A c-replica is a correlated set of PDF replicas, all obtained by determining the best-fit θ (k) Eq. (2.4) but with different (fixed) values of α s : given the data replica D (k) , the minimization Eq. (2.4) is performed several times, for a range of fixed values of α s (m Z ). A c-replica thus corresponds to as many standard NNPDF replicas as the number of values of α s for which the minimization has been performed, all obtained by fitting to the same underlying data replica D (k) .
One can then determine the best-fit value α (k) s for the kth c-replica by minimizing as a function of α s the figure of merit χ 2 Eq. (2.3) computed with θ (k) (α s ) as α s is varied for fixed k. Namely, we first define the figure of merit computed for the kth c-replica, which we can view as a function of α s . Note that the dependence of the theory prediction T and thus of the figure of merit Eq. (2.3) on α s is both explicit, and implicit through the best-fit parameters θ (k) (α s ). We then determine the bestfit value of α s for the kth c-replica as Note that while, as discussed above, in order to avoid overfitting, the best-fit θ (k) is not the absolute minimum of the figure of merit, no overfitting of α s is possible, because overfitting happens when fitting a function, not a single parameter. Hence, in Eq. (2.6) the best fit value α (k) s does denote the absolute minimum. Therefore, in practice α (k) s can be determined by fitting a parabola to the discrete set of values of χ 2 (α s ) for each replica, and finding the minimum of the parabola.
Note also that determining the best-fit for the kth c-replica by first minimizing with respect to θ and then minimizing with respect to α s is equivalent to simultaneously minimizing in the (α s , θ) hyperspace, provided the same figure of merit is used for PDF and α s determination. For instance, the absolute minimum in (α s , θ) is the solution to the coupled equations where Eq. (2.7) is actually a system of N par equations because θ is an N par -component vector and the partial derivative is a gradient. On the other hand, this solution can also be found (compare Fig. 1) by first finding the solution θ(α s ) to Eq. (2.7), determining χ 2 (α s ) = χ 2 (α s , θ(α s )), and solving This two stage procedure yields the same solution as the coupled Eqs. (2.7)-(2.8) because the second term in brackets on the r.h.s. of Eq. (2.9) vanishes since θ(α s ) was the solution of Eq. (2.7). One thus ends up, for each data replica D (k) , with a best fit value (α (k) s , θ (k) ) of both α s and the PDF parameters. That is, from each c-replica we extract a single best fit value α (k) s -an "α s replica" -exactly on the same footing as all the other fit parameters. The ensemble of values α (k) s obtained from all the c-replicas then provides a representation of the probability density of α s from which we can perform statistics in the usual way. Interestingly, this means that we can now not only compute the best fit α s and its uncertainty as the mean and standard deviation (or indeed 68% confidence interval) using the α s replicas, but also the correlation between α s and the PDFs or indeed any PDF-dependent quantity.
In summary, the correlated replica method is akin to the standard NNPDF methodology in that it starts by producing a set of replicas of the original data, but uses these to construct a set of correlated α s -dependent PDF replicas, the c-replicas, which correspond to parameters θ (k) (α s ) when k runs over the replica sample and α s takes a number of discrete values. From each c-replica a best-fit α (k) s can then be determined, so each c-replica yields an α s replica, with α (k) s defined by Eq. (2.6).
Hence, the correlated replica method exploits the fact that in the NNPDF approach it is sufficient to know the best-fit set of parameters for each replica, and all other information is contained in the replica sample. The price to pay for this is that the statistics of the α s fitting is inevitably more demand-ing than with the method of Refs. [17,18] because we have now have to fit a different parabola for each c-replica. The issues arising from this will be discussed in the next section.

Implementation
Building on the conceptual strategy described above, we now present the practical implementation of the correlated replica method. As already mentioned, the best-fit α (k) s Eq. (2.6) for the kth c-replica is determined by fitting a parabola to the figure of merit χ 2 (α s ), viewed as a function of α s , known at the discrete set of α s values for which best-fit θ (k) (α s ) are available. The reliability of the quadratic approximation to χ 2(k) Eq. (2.5) and the stability of the position of the minimum upon inclusion of higher order terms can be studied using standard methods and will be discussed in Sect. 3.2 below.
The best-fit α s and its uncertainty are then determined, according to standard NNPDF methodology, as the mean and standard deviation computed over the sample of α s replicas s is given by Eq. (2.6). The uncertainty due to the finite size of the replica sample can be estimated by bootstrapping. To this purpose, one constructs N res resamples of the original sample of N rep values α (k) s . Each resample is obtained by drawing at random N rep values from the original sample by allowing repetition. This means that each resample differs from the original sample because some values are repeated and others are missing. The finite-size uncertainty is then estimated by first computing the mean α (res,i) s for each of the resamples, where the mean is computed over the N rep values of the ith resample. The bootstrapping estimate of the finite-size uncertainty on the central value of α s is then the standard deviation of the set of α (res,i) The uncertainty on the uncertainty σ can be similarly computed by first determining the uncertainty Eq. (2.10) for each resample, thus leading to an uncertainty σ (res,i) α , and then computing the standard deviation of the ensuing uncertainties: We find that results become independent of the random seed used to generate the bootstrapping resamples when N res 10,000. It turns out that, when determining the best-fit θ (k) (α s ) through the standard NNPDF minimization algorithm, a certain amount of fluctuation of individual values of χ 2 (α s ) about the parabolic best-fit is observed. In other words, the χ 2 profiles as a function α s are not very smooth. It is therefore advantageous to introduce an improvement of the algorithm, called batch minimization, which increases its accuracy at the cost of increasing the time required for fitting.
Furthermore, when using the standard NNPDF minimization, occasionally the fit fails to satisfy a number of convergence and quality criteria (see Sect. 3.3.2 of Ref. [26]), in which case it is discarded. Consequently, for some c-replicas χ 2 (α s ) is not available for all α s values. One must then decide on a sensible criterion for c-replica selection, with the most restrictive criterion being to only accept c-replicas for which all χ 2 (α s ) values are available, and the least restrictive one to accept c-replicas for which at least three χ 2 (α s ) values are available so a parabola can be fitted. We now discuss batch minimization and replica selection criteria in turn.
The idea of batch minimization is to refit a given set of data replicas more than once. In order to improve the smoothness of the χ 2 profiles obtained by the direct use of NNPDF minimization, we exploit the fact that the minimization algorithm is not deterministic, and thus simply rerunning the minimization from a different random seed leads to a slightly different answer. Each of these refits is called a batch. For each c-replica k and each α s value we then end up with several best-fit results θ (k) i (α s ), where i runs over batches. We then pick for each c-replica k and for each α s value the batch which gives the best χ 2 . We also impose the condition that at least two of the batches for the given c-replica and α s value have converged, in order to mitigate the influence of outliers that narrowly pass the post-selection fit criteria. The dependence of results on the number of batches used can then be assessed a posteriori by comparing results found with different numbers of batches.
After batch minimization, we end up with a set of creplicas θ (k) (α s ) where, however, for several c-replicas, results may be missing for one or more α s values because convergence was not achieved. We must thus determine the minimum number of α s values N min such that a c-replica is accepted. The threshold N min is chosen to ensure the stability of results. Curves with too few points lead to an unreliable parabolic fit, and thus an unreliable best-fit α (k) s for that c-replica. This then leads to outlier values of α (k) s and a spuriously large value of the uncertainty on the α (k) s determination. On the other hand, once the number of points is sufficient for a reliable parabolic fit, requiring more points does not improve the determination of α (k) s , but it reduces the number of c-replicas which are retained in the final sample, which in turn increases the finite-size uncertainty.
Therefore, the optimal value of N min arises from a tradeoff between the uncertainty on α (k) s from the parabolic fitting, and the finite-size uncertainty. In order to keep both criteria into account, we fix N min by minimizing the bootstrapping uncertainty σ Eq. (2.13). However, in order to make sure that the selection is not too tight, we do not minimize σ itself. Rather, we first multiply it by a penalty factor that depends on the number of points. This is in turn determined as the 99% confidence level factor from a two sided Student t distribution. Indeed, if the distribution of best-fit α (k) s over replicas is Gaussian, then the difference between the sampled and true central value follows a Student t distribution with N rep − 1 degrees of freedom, zero mean and scale parameter σ / N rep . A given confidence level around the mean is equal to the standard deviation σ T CL,(N rep −1) , where T CL,N is the percentile at CL confidence level for the two-sided confidence factor of the Student t distribution with N degrees of freedom. Hence, we choose a 99% confidence level, and we determine N min as the value yielding the minimum of σ T 0.99,(N rep −1) . Also in this case, the dependence of results on the choice of selection criteria can be studied a posteriori, and will be discussed in Sect. 3.2.

The strong coupling constant from NNPDF3.1
We now present the main result of this work, namely the determination of α s (m Z ) based on the methodology discussed in Sect. 2. We first present the best-fit result for α s and its experimental uncertainty, determined through the correlated replica method. We then discuss methodological and theoretical uncertainties. We finally collect our final result and briefly compare it to other recent determinations from PDF fits and to the PDG average.

Best-fit results for α s and statistical uncertainty
We have determined α s (m Z ) both at NLO and NNLO using the methodology and dataset of the NNPDF3.1 global analysis [16]. The only difference in the fit settings is the theoretical description of the inclusive jet production datasets at NNLO. Here we use exact NNLO theory [33] for the ATLAS [34] and CMS [35] inclusive jet measurements at 7 TeV, and discard the other jet datasets used in NNPDF3.1 for which the NNLO calculation is not available (note that, as in NNPDF3.1, only ATLAS data in the central rapidity bin are included). To ensure a consistent comparison, the input datasets of the NLO and NNLO fits are identical, up to small differences in the kinematical cuts as explained in [16].
Specifically, we determine α s by generating a set of 400 data replicas, and from them a set of 400 c-replicas each with 21 values of α s , thus corresponding to a total of 8400 PDF replicas correlated as discussed in Sect. 2.1. These creplicas are generated for α s (m Z ) ranging between 0.106 and 0.130, varied in steps of α s = 0.002 between 0.106 and 0.112 and between 0.128 and 0.130, and in steps of α s = 0.001 between 0.112 and 0.128, adding up to the total of 21 values. From these we determine α s replicas, which form a representation of the probability distribution of α s .
At NNLO we find This result is based on a total of N rep = 379 c-replicas, selected from a starting set of 400 after batch minimization of three batches, using the minimization and selection methods described in Sect. 2.2. At NLO we find In this case, the sample includes N rep = 108 c-replicas selected after batch minimization with two batches. The smaller number of c-replicas selected at NLO is in part explained by the requirement (see Sect. 2.2) that two batches have converged for the given α s value, which is of course less severe when three batches are available, but the worse quality of the NLO fit also plays a role since it causes more fits to be discarded by the post-selection criteria. The uncertainty quoted in Eqs. (3.1) and (3.2) is that obtained using standard NNPDF methodology, namely, taking the standard deviation over the α s replica sample. We have verified that essentially the same results are obtained if instead we compute the 68% confidence interval. The uncertainty is obtained in precisely the same way as our PDF uncertainty, to which it is strongly correlated; it includes the propagated correlated uncertainty from the underlying data, and uncertainties coming from possible inefficiencies of the minimization procedure. This uncertainty is what we refer to as the experimental uncertainty on α s (m Z ). It will have to be supplemented by methodological and theoretical uncertainties, to be discussed in Sects. 3.2 and 3.3 below.
The 379 c-replicas selected for the NNLO determination are shown in Fig. 2. The color scale of each curve indicates the best-fit α s value. It is apparent that the vast majority of the curves exhibit an approximately parabolic behaviour. The probability distributions of the best-fit values α (k) s Eq. (2.6) which correspond to each c-replica, both at NLO and at NNLO, are shown in Fig. 3, where the markers indicate the value of α (k) s for each specific c-replica. These probability densities have been determined using the Kernel Density Estimate method, see [36]. We find that the probability distribution for α s (m Z ) is both shifted to higher values and broadened when going from NNLO to NLO. The decrease of the best-fit value of α s (m Z ) when going from NLO to Fig. 2 The χ 2 profiles for each of the 379 c-replicas used for the NNLO determination of α s (m Z ), Eq. (3.1). Each curve corresponds to an individual c-replica, and the color scale indicates the best-fit α s value determined from the parabolic fit to that curve NNLO has been repeatedly observed before (see Table 1 of Ref. [37] for an extensive set of examples), also in our previous determination [17,18], while the broadening is due to the poorer quality of the NLO fit.
The impact on the α s determination of any subset of the input data can be roughly assessed by studying its contribution to the total figure of merit. We have done this by determining replica by replica the corresponding partial χ 2 p for a process (or group of processes) p, defined as the figure of merit Eq. (2.3) with the summation over i, j now restricted to data which belong to the specific subset p. The α s fit procedure through the correlated replica method is then just repeated but using this partial χ 2 p . Namely, for each c-replica the partial χ 2(k) p for process p is computed, a parabola is fitted to it, the corresponding minimum α (k) s, p of the parabola is determined, and the resulting set of minima is used to find the value of α s (m Z ) and its uncertainty. Here we consider the following eight groups of processes p: top production, the Z p T distributions, collider and fixed target Drell-Yan, inclusive jets, and deep-inelastic scattering (DIS) either at HERA or at fixed-target experiments, in the latter case separating charged lepton and neutrino beams. The number of data points corresponding to each of these data subsets is shown in Table 1. Not unexpectedly, the χ 2(k) p profiles for data subsets turn out to be rather less parabolic than the total χ 2 , especially for processes such as neutrino DIS or fixed target Drell-Yan that have weak sensitivity to α s .
When determining α s (m Z ) from the partial χ 2(k) p , we do not repeat the replica selection and simply use the same replicas selected for the total dataset. Consequently, we must apply a form of post-selection, whereby each time a parabola for χ 2(k) p has no minimum the corresponding c-replica is ignored. At NNLO, for five out of eight data subsets we retain all 379 c-replicas, while for jets, neutrino DIS, and fixed-target Drell-Yan, we retain only 376, 366, and 302 c-replicas respectively. The results for the partial α s (m Z ) determined from χ 2 p for the various families of processes are   Table 1 at NLO and NNLO collected in Fig. 4. The central value and uncertainty shown are respectively determined as the median and 68% symmetric confidence level interval from the corresponding partial α (k) s, p . This is because the analogue of Fig. 3 for individual processes turns out to be rather non-gaussian, especially for processes such as fixed-target Drell-Yan that only have a weak handle on α s .
The values of α s (m Z ) shown in Fig. 4 should be interpreted with some care. Indeed, the partial χ 2 p is in each case computed using PDF c-replicas determined from the minimization of the global χ 2 . These are in general different from the c-replicas that would be determined by simultaneous minimization of χ 2 p with respect to α s and the PDFs. Therefore, the values of α s, p in Fig. 4 cannot be interpreted as the best-fit values of α s (m Z ) for a given subset p. They instead provide an estimate of the pull on the best-fit α s (m Z ) value that specific families of processes have within the global fit subject to the constraints from the rest of the data.
Moreover, even their interpretation as pulls is only approximate. Firstly, the replica selection is applied to the total χ 2 rather than to each partial χ 2 p , so that several partial χ 2(k) p profiles turn out not to have a minimum. Furthermore, the total χ 2 includes cross-correlations which are lost when determining partial χ 2 p , because the covariance matrix C t 0 in Eq. (2.3) is generally nonzero even when i and j belong to different data subsets. For instance, inclusive jet, Z p T , and Drell-Yan measurements from the same experiment (ATLAS, or CMS) are correlated amongst themselves by the common luminosity uncertainty. Finally, partial α s values are correlated through the underlying PDFs, implying that the pulls should not be expected to combine additively into the final result.
Even with all these caveats, Fig. 4 shows that the very accurate α s (m Z ) value from the global dataset is obtained from a combination of pulls which correspond to values of α s (m Z ) dispersed about the global best-fit value, without signs of tension or inconsistency, and subject to significant fluctuations which are suppressed when constructing the total χ 2 . This supports our conclusion that the current determination of α s (m Z ) from a global fit is more precise and accurate than determinations based on subsets of data relying on preexisting PDF sets.
Finally, we compare the current NNLO determination of α s (m Z ), Eq. (3.1) and Fig. 4, with the one found using the method of Refs. [17,18]. We fix α s and add the contribution to the χ 2 from each standard PDF replica for that α s value. We then determine the total χ 2 (α s ), fit a parabola to it, and determine the best-fit and uncertainty as the minimum and χ 2 = 1 interval. For simplicity, we do this without using batch minimization, i.e. we compute the total χ 2 from one of the batches (batch II, see Sect. 3.2 below) which then enter the batch minimization procedure. Using this method we find Also in this case we can repeat the determination for different data subsets based on the partial χ 2 p , and the corresponding results are compared in Fig. 5.
As expected, and discussed in the introduction and in Sect. 2.1, we find that the best-fit values of α s (m Z ) determined with the old method [17,18] and with the new correlated replica method are in good agreement, both for the global dataset and for the data subsets. The small differences in central values are most likely due to uncertainties related to the finite size of the replica sample, which, as discussed in [17,18], can be non-negligible when the old method is used. On the other hand, also as expected, Fig. 5 Comparison of the NNLO determination of α s (m Z ) using the method of [17,18], which neglects the correlation between α s and PDFs, and the current one based on the correlated replicas Fig. 6 The NNLO cumulative differences, χ 2 p (α s ) − χ 2 p (0.1185), between the partial χ 2 p values evaluated at α s (m Z ) and at best-fit value α s (m Z ) = 0.1185 for different families of processes. In the part of the plot above (below) zero, only contributions from experiments for which the cumulative difference is positive (negative) are shown (see text). The plot is displayed either with a wider (left) or narrower (right) choice of range on the y axis neglecting the correlation between α s and PDFs as in the old method leads in general to an underestimate of the uncertainty on α s . This effect is more marked for processes such as fixed-target Drell-Yan and neutrino DIS that have a limited sensitivity to α s , because in this case the difference in length of the semi-axes of the error ellipse in Fig. 1 is large.
This determination of α s (m Z ) from the total χ 2 also offers a complementary way of quantifying how much each family of processes constrains the final best-fit value, by plotting the contribution of each data subset to the total χ 2 . Specifically, we show in Fig. 6 the cumulative differences at NNLO, χ 2 p (α s ) − χ 2 p (0.1185), between each partial χ 2 p and its value computed at the global best-fit α s (m Z ) value, neglecting cross-correlations between different data subsets. The plot is divided into two halfs: above zero, only positive differences are shown, and below zero, only negative ones. Thus, when all differences are positive the plot shows the breakdown of the total χ 2 into the contribution of different experiments (up to neglected cross-correlations), while when some of them are negative the lower part of the plot shows by how much the χ 2 of the individual experiments shown has improved in comparison to their value at the global minimum α s (M z ) = 0.1185). In order to increase readability, the plot is displayed twice, with two different choices of scale on the y axis.
From this comparison, we observe that the LHC data significantly contribute to constraining α s . In particular, it is interesting to note that the 13 data points from topquark pair production lead to a significant contribution to the total χ 2 away from the best-fit, even though the global dataset contains almost 4000 data points. Similar considerations apply to the Z p T distributions. This means that there is a small range of values of α s where these two groups of processes are consistent with the rest of the data entering the fit, thereby providing a tight constraint on α s .

Methodological uncertainties
In view of the rather small experimental uncertainty on the final value of α s (m Z ), Eqs. (3.1)-(3.2), we need to assess possible uncertainties associated to the various aspects of our methodology described in Sect. 2. Specifically, we discuss here the methodological uncertainties associated to c-replica selection, batch minimization, the quadratic approximation to χ 2 profiles, and the treatment of correlated systematics.
The replica selection algorithm determines an optimal value of N min , the minimal number of α s for which results must be available for a given c-replica to be selected. We have varied this value from its minimum N min = 3 (needed in order to fit a parabola) to a high value N min = 18 (meaning that at most three values α s can be missing in order for a c-replica to be retained). Results for the number of c-replicas passing the criterion and the ensuing value of α s are collected in Table 2 for a number of choices. In each case we also show the finite-size uncertainty α s on the best-fit α s estimated by bootstrapping, Eq. (2.12).
The number of surviving c-replicas varies significantly; all the starting 400 c-replicas pass the loosest criterion (i.e., it is always possible to fit a parabola to any c-replica), but only N rep = 12 c-replicas pass the most restrictive criterion. However, even with this most restrictive criterion the finite-size uncertainty is below the permille level. For the value selected by the algorithm, the finite-size uncertainty is of order 0.03%, i.e. by almost a factor 20 smaller than the experimental uncertainty Eq. (3.1) and it does not decrease further even when all c-replicas are kept. The finite-size uncertainty on the α s uncertainty σ itself Eq. (2.13) is comparable in all cases.
The value of α s (m Z ) and its experimental uncertainty are hence very stable; the shift of central value and uncertainty when the selection criterion is varied is always smaller than the finite-size uncertainty. This stability can be understood by observing that each c-replica consists of at least N min correlated PDF replicas, so each of the determinations shown in Table 2 is obtained from more than N min × N rep PDF replicas. We thus estimate that the bootstrapping uncertainty, and the related but smaller uncertainty due to choice of replica selection, to be of order α s = 0.00003 (0.03%), one order of magnitude smaller than the experimental uncertainty.
We next turn to discuss batch minimization. The results shown in Table 2 all correspond to the NNLO baseline which uses batch minimization with three batches. In order to assess the impact of batch minimization, in Table 3 we compare results obtained with each of the three batches, with the three possible pairs, and combining the three batches. In each case we show the final best-fit α s (m Z ) and experimental uncertainty, the value of N min , the minimum number of α (k) s values per c-replica, and the number of surviving c-replicas N rep .
It is clear from this comparison that as more batches are combined, results become more stable. The values of N min are on average larger with two batches, and larger still with three, but without a reduction of the number of surviving creplicas N rep as was observed in Table 2. With three batches, N rep is largest even though N min is also largest. This means that, thanks to batch minimization, the number of available α (k) s values per replica is on average higher. It follows that the finite-size uncertainty is reduced by batch minimization, thus leading to the very small uncertainties shown in Table 2.
The values of α s (m Z ) behave as expected upon use of batch minimization. The experimental uncertainty is reduced when more batches are used and the central values with different combinations of batches are all consistent with each other within given uncertainties. Furthermore, the differences in central values with different combinations of batches are reduced upon use of batch minimization (they are smaller when using two batches than when using a single batch). Table 4 Results for the NNLO determinations of α s (m Z ) when the N trim outer values of α s are not used and the fit is restricted to a smaller range. In the bottom part of the table we also show results found discarding values asymmetrically, at the upper or lower edge of the range. In each case we show the number of discarded α s values, the best-fit value of α s (m Z ), and the number of surviving c-replicas N rep . The first row (in boldface) corresponds to our final result Eq. (3.1) Additionally, the shift in central value when increasing the number of batches is rather smaller than the uncertainty, and, finally, the central value is stabilized when increasing the number of batches, so the difference between two and three batches is on average smaller than the difference between one and two batches. We conclude that the value of α s (m Z ) found using three batches is the most accurate. We observe that even the shift between the three-batch value and the single-batch value which differs most from it is about a third of the finite-size uncertainty. We take this as further evidence that there is no extra contribution of methodological origin due to batch minimization to be added to the statistical uncertainty. We finally observe that the two-batch result is in fact consistent within its very slightly larger uncertainty, thus justifying the use of only two batches at NLO.
We next turn to the methodological uncertainties related to the quadratic fitting of χ 2 profiles. We have studied this in three different ways: by removing outer values of α s (m Z ) from the fit; by adding higher order terms to the fitting function; and by changing the fitting variable. We discuss each in turn.
First, we have repeated the NNLO determination removing α s values that are farthest from the best-fit value α s (m Z ) = 0.1185, fitting a smaller range of values around the minimum. As a further consistency check, we have removed α s values asymmetrically. Results are shown in Table 4; in each case we show the number of discarded α s values N trim , the resulting fitted range, the best fit α s (m Z ) and uncertainty, and the number of surviving c-replicas N rep . Here too, the behaviour is consistent with expectations. As the fitted range is reduced, the experimental uncertainty increases and the number of surviving c-replicas decreases (thereby also increasing the finite-size uncertainty). The central value, however, is extremely stable; the shift in central Table 5 Same as Table 2, comparing the default parabolic fitting (in boldface) of the χ 2 (α s ) profiles with those with a transformed input, both χ 2 (ln(1 + α s )) and χ 2 (exp(α s )) value when restricting the range is always more than a factor two smaller than the experimental uncertainty. In fact, the shift is never larger than = 0.00010 (0.08%) unless the number of surviving c-replicas becomes of order ten, in which case the finite-size uncertainty (recall Table 2) is of the same order or larger.
A different way of testing for deviations from quadratic behaviour is to apply a criterion to assess fit quality to both quadratic and cubic fits. Here we use the Akaike Information Criterion (AIC) [38], which estimates the expected relative distance between a given fitted model and the unknown underlying law [39]. The AIC score balances goodness of fit against simplicity of the model. A lower score corresponds to a lower expected distance measured by the Kullback-Leibler divergence. The AIC score is defined by where r is the number of degrees of freedom of the model, n is the number of fitted points, and ln(L) is the log-likelihood associated with the model. In our case, we fit to χ 2(k) (α s ), Eq. (2.5), viewed as a function of α s using either a parabola (as in our default determination) or a higher order polynomial. The log-likelihood is then in each case just the χ 2 of this fit. Computing the AIC score for each fitted profile, averaging over c-replicas, and taking the variance of results as a measure of the uncertainty, we find AIC = 169 ± 37 for the default quadratic fit and AIC = 173 ± 35 for a cubic fit. We conclude that there is no evidence that a cubic fit is better than a quadratic one.
We perform a final test based on the observation that any transformation of the error function profile of the form (3.5) where f is sufficiently smooth and monotonic, should lead to the same best-fit value of α s . The results of fitting α s from the transformed profiles Eq. (3.5) with f (α s ) = exp(α s ) and f (α s ) = ln(1 + α s ) are shown in Table 5. The argument of the log is shifted so that f (α s ) admits a Taylor expansion in powers of α s . Reassuringly, we find extreme stability with respect to these transformations of the fitting argument. Combining results from Tables 4 and 5 and the analysis based on the AIC score we can conservatively take as an estimate of the uncertainty related to parabolic fitting the largest shift observed in Table 2, neglecting the cases with N rep < 100 which are dominated by finite-size uncertainty, namely par = 0.00010 (0.08%). (3.6) We finally turn to the uncertainty related to the treatment of experimental correlated systematic errors. As mentioned in Sec. 2.1, the covariance matrix in the presence of multiplicative uncertainties should not be identified with the experimental covariance matrix, in order to avoid biasing the fit [23]. We thus adopt the t 0 method, introduced in [24], benchmarked in [25], and used for the determination of all NNPDF sets from NNPDF2.0 [40] onwards. In this procedure, the normalization of the multiplicative uncertainties that enter the covariance matrix is iteratively determined from a prior the-ory prediction. Because the PDFs and α s are now determined on the same footing, the same covariance matrix is used for both. It is clear that the same χ 2 definition must be used in Eq. (2.9) as in Eqs. (2.7)-(2.8) in order for the same minimum to be found.
Indeed, it is interesting to note that using an inconsistent definition of the covariance matrix significantly biases the result of α s (m Z ). In Fig. 7 we compare the distribution of NNLO α s (m Z ) values as well as the total and partial best-fit values and uncertainties, computed for a single batch, either consistently using the t 0 covariance matrix (see Figs. 3, 5 for the corresponding results with three batches) or inconsistently using the experimental covariance matrix. We find that the inconsistent definition leads to a much broader distribution for the total χ 2 , thereby signaling the lack of consistency, and, more importantly, a biased central value α s (m Z ) = 0.114 ± 0.001 exp (0.9%), shifted by about 9σ in comparison to the correct result Eq. (3.1). The fact that a downward shift of α s (m Z ) is observed when using the inconsistent definition can be understood based on the observation that the bias [23] typically leads to the best-fit undershooting the data, essentially because with multiplicative uncertainties a lower prediction has a smaller uncertainty [41]. Indeed, inspection of the partial best-fit values shows that the bias is much stronger for collider experiments than the fixed-target ones. This is what one would expect, because systematic uncertainties are multiplicative for collider experiments, while they are mostly additive for fixedtarget [25], so any effect or bias related to the treatment of multiplicative uncertainties should be mostly seen in collider data. The use of the t 0 procedure in principle leads to a further methodological uncertainty related to the choice of the prior used for the construction of the t 0 matrix, which should therefore be assessed. In order to determine the final result Eq. (3.1) the t 0 matrix was constructed using the best-fit PDF set from batch II of Table 3. We have repeated the determination constructing the t 0 matrix from the best-fit PDF set of either of the other two batches. Results are collected in Table 6. It is clear that, using the consistent t 0 method, results are extremely stable. We can conservatively estimate the uncertainty due to the choice of t 0 from the largest shift seen in Table 6 as t 0 = 0.00004 (0.03%).
In summary, we conservatively estimate methodological uncertainties by adding in quadrature the finite-size uncertainty α s = 0.00003, the uncertainty related to the parabolic approximation par = 0.00010 and the uncertainty related to the treatment of correlated systematics t 0 = 0.00004, with the result σ meth = 0.00011 (0.09%). (3.7) Therefore, we find that, at NNLO, methodological uncertainties are smaller than the experimental uncertainties Eq. (3.1) by a factor five.

Theoretical uncertainties from missing higher orders
A determination of α s (m Z ) is dependent on the perturbative order of the QCD calculations on which it relies. Therefore, at any fixed order it is affected by a missing higher order uncertainty (MHOU). In older, and also some more recent determinations of α s (m Z ) (specifically for determination in PDF fits see Refs. [17,42,43]) no attempt was made to estimate the MHOU, and sometimes NLO or NNLO values of α s (m Z ) were quoted with the understanding that they might differ by an amount greater than the quoted uncertainty due to this missing uncertainty. However, as the experimental uncertainty decreases, an estimate of the MHOU becomes mandatory, and in the context of PDF fits it was done e.g. in Ref. [18]. Indeed, this uncertainty, usually estimated by scale variation, is typically dominant in more recent determinations [9][10][11][12][13][14][15].
In the present case, a first handle on the MHOU associated to α s is provided by the difference between the NLO and NNLO results Eqs. which corresponds to a 2% shift of the NNLO central value. This is about four times larger than the experimental uncertainty in Eq. (3.1), thereby suggesting that even at NNLO the MHOU on the α s (m Z ) determination might be comparable to, or larger than the experimental uncertainty.
In our previous determination of α s Ref. [18] the MHOU was estimated using the Cacciari-Houdeau (CH) method [44], which relies on a Bayesian estimate of the missing higher perturbative orders based on the behaviour of the known orders. Use of exactly the same method of Ref. [18], to which the reader is referred for details, leads to the values CH, NLO = 0.003, (3.9) CH, NNLO = 0.0004, (3.10) for the 68% confidence level MHOU on α s (M Z ). The rather large difference in the MHOU estimate between NLO and NNLO stems from the fact that there is a significant shift when going from LO to NLO, but a much smaller one when going from NLO to NNLO. The NLO estimate of the MHOUs in Eq. (3.9) is reassuringly in good agreement with the observed shift Eq. (3.8). The NNLO uncertainty Eq. (3.10) is also consistent with expectations based on the CH uncertainty estimate of Ref. [18], where the value of α s (m Z ) determined using the NNPDF2.1 set was found to lead to CH, NNLO = 0.0009. Indeed, PDF uncertainties in the NNPDF3.1 set are generally smaller than those on NNPDF2.1 by a factor of two or more, due to significant impact of LHC data in the more recent determination.
In addition, the shift between NLO and NNLO PDFs is found to be smaller in NNPDF3.1 than in previous NNPDF sets [45], presumably because MHO terms pull in different directions and thus partly cancel each other to a greater extent in a more global fit. Indeed, we find a similar increase of perturbative stability of PDFs and of the associated α s (m Z ) by repeating the analysis presented here for reduced datasets [46]. Therefore, the reduction of the MHOU by a comparable factor in Eq. (3.8) in comparison to Ref. [18] is expected.
Nevertheless, the very small value of the MHOU at NNLO, Eq. (3.10), even smaller than the already small experimental uncertainty Eq. (3.1), may seem rather too optimistic. There are furthermore several reasons of principle and practice why the reliability of the CH method in the present case is dubious. The main one is that the implementation of the method suggested in Ref. [18] relies on a guess for an underlying "true" value α We therefore prefer to adopt a more conservative estimate. Namely, we assume that the MHOU on the NNLO result is half the difference between the NLO and NNLO results Eq. (3.8): about twice the size of the corresponding experimental uncertainty Eq. (3.1). Whereas this is surely a very crude estimate, we do not feel that any of the available methods can lead to a more reliable conclusion. On top of the missing higher fixed-order QCD corrections, several other aspects of the theory used in the simultaneous determination of α s (m Z ) and PDFs also lead to uncertainties. These include the values of the heavy quark masses, standard model parameters (specifically CKM matrix elements and electroweak couplings), electroweak corrections, QCD resummation corrections [47,48], QCD power corrections, and nuclear corrections. Many of these uncertainties were assessed in the NNPDF3.1 PDF determination that we are relying upon [16], and found to be smaller than PDF uncertainties. In particular, the dependence on the charm mass in previous PDF determinations is substantially reduced in NNPDF3.1 and likely rather smaller than the MHOU, thanks to the presence of an independently parametrized charm PDF [49], and electroweak corrections are carefully kept under control thanks to the choice of suitable kinematic cuts. But PDF uncertainties mix with the experimental uncertainty on α s (m Z ), with which they are strongly correlated, and are in fact indistinguishable from it, as discussed in Sect. 2.1, so the hierarchy of uncertainties on PDFs and α s (m Z ) is the same. We conclude that we have evidence that most of these theoretical uncertainties are subdominant in comparison to the experimental uncertainty Eq. (3.1), and thus even more so in comparison to the MHOU Eq. (3.11).

Final results and comparisons
We can now collect results. Combining the NNLO value and experimental uncertainty Eq. (3.1), the methodological uncertainty Eq. (3.7) and the theoretical uncertainty Eq. (3.11) we get α NNLO s (m Z ) = 0.1185 ± 0.0005 exp ± 0.0001 meth ± 0.0011 th = 0.1185 ± 0.0012 (1%), (3.12) where in the last step we have added all uncertainties in quadrature. For a comparison to other determinations, such as the PDG average, we recommend using only the experimental uncertainty (the methodological uncertainty being negligible), which reflects the limitations of our result and procedure, but not the limitation due to the fact that our result is obtained at NNLO. For precision phenomenology, however, we recommend use of the total uncertainty in order to conservatively account for the MHOU.
This result can be compared to the previous one [18] based on NNPDF2.1, α NNLO s (m Z ) = 0.1173 ± 0.0007 exp ± 0.0009 th . In comparison to this older result, the central value of α s (m Z ) has increased by α s = +0.0012 . As far as uncertainties are concerned, both the theoretical and experimental uncertainties on this previous result are larger, if one compares like with like. The experimental uncertainty should actually be compared to Eq. (3.3) as it was obtained with the same method. The uncertainty is somewhat underestimated because it neglects the correlation between PDFs and α s , while the theory uncertainty should be compared to Eq. (3.10) which is also based on the CH method.
We conclude that, in comparison to Ref. [18], the current result is more precise, though with more conservatively estimated uncertainties.
In Fig. 8 we compare the NNLO result of Eq. (3.12) to our previous result [18], to the current PDG average [3], and to two recent determinations obtained from simultaneous fit of PDFs and α s (m Z ), ABMP16 [43] and MMHT2014 [42]. We find good agreement with the PDG average as well as with the MMHT14 and NNPDF2.1 determinations. It has been suggested [50,51] that the lower ABMP16 value can be partly explained by the use of a fixed-flavour number scheme with N f = 3 for the treatment of DIS data. It is interesting to observe that the current AMBP16 value is higher than previous values of α s (m Z ) obtained by the same group [52], from which the ABMP16 analysis in particular differs because of inclusion in Ref. [43] of LHC top production and W and Z production data (described with N f = 5).
Interestingly, the α s (m Z ) determination from the NNPDF3.1 fit is higher than any other recent determination from PDF fits. Inspection of Figs. 4 and 6 strongly suggests that this increase is driven by the high-precision LHC data, especially for gauge boson production (including the Z p T distribution) but also for top and jet production.

Summary and outlook
In this work we have presented a new determination of the strong coupling constant α s (m Z ) jointly with a global determination of PDFs which, by relying on NNPDF3.1, for the first time includes a large amount of LHC data using exact NNLO theory in all cases. In comparison to a previous determination based on NNPDF2.1, our results exploit the new correlated replica method that is equivalent to the simultaneous fit of PDFs and α s . This new method thus fully accounts for the correlations between PDFs and α s in the determination of the best-fit value of α s and of the associated uncertainty.
We find that the determination of α s (m Z ) is considerably stabilized by the use of a wide set of different processes and data, and we provide evidence that a global simultaneous determination of α s (m Z ) and PDFs leads to a more stable and accurate result than the one obtained from subsets of data. We thus obtain a value of α s (m Z ) which is likely to be more precise and more accurate than previous results based on similar techniques. We find that the LHC data consistently lead to an increase in the central value of α s (m Z ), and observe good overall consistency between the datasets entering the global fit. Our NNLO determination turns out to be in agreement within uncertainties with previous results from global fits and with the PDG average.
The main limitation of our result comes from the lack of a reliable method to estimate the uncertainties related to missing higher order perturbative corrections. Theoretical progress in this direction is needed, and perhaps expected, and would be a major source of future improvement. For the time being, even with a very conservative estimate of the theoretical uncertainty, our result provides one of the most accurate determinations of α s (m Z ) available, and thus provides valuable input for precision tests of the Standard Model and for searches for new physics beyond it.