Determination of the strong coupling constant (cid:11) S ( m Z ) from measurements of inclusive W (cid:6) and Z boson production cross sections in proton-proton collisions at p s = 7 and 8 TeV

: Twelve measurements of inclusive cross sections of W (cid:6) and Z boson production, performed in proton-proton collisions at centre-of-mass energies of 7 and 8 TeV, are compared with perturbative quantum chromodynamics calculations at next-to-next-to-leading order (NNLO) accuracy obtained with the CT14, HERAPDF2.0, MMHT14, and NNPDF3.0 parton distribution functions (PDFs). Data and theory agree well for all PDF sets, taking into account the experimental and theoretical uncertainties. A novel procedure is employed to extract the strong coupling constant at the Z pole mass from a detailed comparison of all the experimental (cid:12)ducial cross sections to the corresponding NNLO theoretical predictions, yielding (cid:11) S ( m Z ) = 0 : 1163 +0 : 0024 (cid:0) 0 : 0031 (CT14), 0 : 1072 +0 : 0043 (cid:0) 0 : 0040 (HERAPDF2.0), 0 : 1186 (cid:6) 0 : 0025 (MMHT14), and 0 : 1147 (cid:6) 0 : 0023 (NNPDF3.0). Using the results obtained with the CT14 and MMHT14 PDFs, which yield the most robust and stable (cid:11) S ( m Z ) extractions, a value (cid:11) S ( m Z ) = 0 : 1175 +0 : 0025 (cid:0) 0 : 0028 is determined.


Introduction
In the chiral limit of zero quark masses, the α S coupling is the only free parameter of quantum chromodynamics (QCD), the theory of the strong interaction between quarks and gluons. Because of its logarithmic decrease with energy (asymptotic freedom), α S is commonly given at a reference scale, often taken at the Z pole mass. Its current value, α S (m Z ) = 0.1181 ± 0.0011, is known with a ±0.9% uncertainty, making it the least precisely known of all interaction couplings in nature [1]. The precision of the strong coupling value plays an important role in all theoretical calculations of perturbative QCD (pQCD) processes involving partons, and currently leads to 3-7% uncertainties in key Higgs boson processes, such as the cross sections for gluon-gluon fusion (gg → H) and associated production with a top quark pair (ttH), as well as the H → bb, cc, gg partial decay widths [2]. As one of the fundamental parameters of the standard model (SM), the uncertainties of the QCD coupling value also dominate the propagated parametric uncertainties in the theoretical calculations of the top quark mass [3], as well as of electroweak (EW) precision observables [4]. Last but not least, α S also impacts physics approaching the Planck scale, -1 -
either through the EW vacuum stability [5], or in searches of new coloured sectors that may modify its running towards the grand unification scale [6,7]. The current α S (m Z ) world-average value is derived from a combination of six subclasses of (mostly) independent observables measured at various energy scales, which are compared with pQCD calculations at next-to-next-to-leading order (NNLO), or beyond, accuracy [1]. The only hadron collider observable so far that provides a constraint on α S at this level of theoretical accuracy is the total tt cross section [8][9][10]. One of the paths towards improvement of our knowledge of the QCD coupling is the inclusion into the world average of new independent observables sensitive to α S that are experimentally and theoretically known with high precision [11,12]. Charged-and neutral-current Drell-Yan processes in their leptonic decay modes, pp → W ± → ± ν and pp → Z → + − with ± = e ± , µ ± , are the most accurately known processes currently accessible in proton-proton (pp) collisions at the CERN LHC. Experimentally, the uncertainties in the inclusive V = W ± , Z production cross sections measured by the CMS experiment are between 3 and 5%; these are dominated by the integrated luminosity uncertainty, whereas the statistical uncertainties are at the subpercent level (table 1) [13,14]. On the theoretical side, the corresponding cross sections are known at NNLO pQCD accuracy [15], with about 1-4% parton distribution function (PDF), and 0.3-1.3% scale uncertainties [16]. Electroweak corrections, which lead to a few percent reduction of the pure-pQCD W ± and Z boson production cross sections, are known at next-to-leading order (NLO) accuracy [17].
Theoretical calculations [18] indicate that about one fourth of the total V production cross sections at LHC energies come from partonic processes beyond the Born level, and thereby depend on the QCD coupling value. By calculating the V production cross sections at NNLO for varying α S (m Z ) values, and by comparing the theoretical predictions to experimental data, one can therefore derive a value of the strong coupling constant at the Z pole, independent of other current extractions [19]. By combining such a result with those derived from other methods, the overall uncertainty in the α S (m Z ) world average -2 -

JHEP06(2020)018
can eventually be reduced. The use of inclusive W ± , Z boson cross sections to extract the QCD coupling is presented here for the first time. This method is similar to the one used to extract α S from the inclusive tt cross sections at hadron colliders [8][9][10], except that the underlying physical process is quite different. Whereas σ(tt) depends on α S already at leading order (LO), albeit with ≈5% theoretical and experimental uncertainties, σ(V) is more precisely known experimentally and theoretically, although at the Born level its underlying partonic processes are purely EW with a dependence on α S that comes only through higher-order pQCD corrections (at LO, the σ(V) cross sections also depend on α S via the PDFs).
The paper is organised as follows. The experimental setup used in the twelve CMS original measurements is summarised in section 2. In section 3, the theoretical tools used to perform the calculations are outlined. In section 4, the experimental and theoretical cross sections with associated uncertainties, are compared. In section 5, the method to extract α S (m Z ) from the data-to-theory comparison for each measurement is described, as well as the approach to combine all α S (m Z ) estimates into a single value per PDF set that properly takes into account the experimental and theoretical uncertainties and their correlations. The final α S (m Z ) values derived are presented and discussed in section 6. The work is summarised in section 7.

The CMS detector
The results presented here are based on a phenomenological study of W ± and Z boson fiducial cross sections measured by the CMS experiment in pp collisions at centre-of-mass (c.m.) energies of √ s = 7 and 8 TeV with integrated luminosities of 38.0 and 18.2 pb −1 , respectively [13,14]. The experimental and theoretical EW boson production cross sections quoted in the whole paper are to be understood as multiplied by their associated leptonic branching fractions, but for simplicity are referred to as "cross sections" hereafter. The final states of interest are those with decay charged leptons (electrons or muons) passing the acceptance criteria listed in table 1.
The central feature of the CMS apparatus is a superconducting solenoid, of 6 m internal diameter, providing a magnetic field of 3.8 T. Within the field volume are a silicon pixel and strip tracker, a crystal electromagnetic calorimeter (ECAL), and a brass and scintillator hadron calorimeter. Electrons with p T > 25 GeV are identified as clusters of energy deposits in the ECAL matched to tracks measured with the silicon tracker. The ECAL fiducial region is defined by |η| < 1.44 (barrel) or 1.57 < |η| < 2.5 (endcap), where η is the pseudorapidity of the energy cluster. Muons are measured in gas-ionisation detectors embedded in the steel flux-return yoke of the magnet. Muons with p T > 20 or 25 GeV and |η| < 2.1 are selected in the analyses. Details of the CMS detector and its performance can be found elsewhere [20].

Theoretical calculations
According to the pQCD factorisation theorem [21], the cross section for the production of a heavy elementary particle in pp collisions can be calculated through the convolution of -3 -JHEP06(2020)018 matrix elements for the relevant parton-parton subprocesses, computed at a given order in the α S expansion evaluated at a renormalisation scale µ r , and a universal nonperturbative part describing the parton density at the factorisation energy scale µ f and parton fractional momentum x i in the proton. The production cross section of an EW boson can be written where the functions f i represent the PDFs of each proton, determined from experimental data, and the expression in brackets is the perturbative expansion of the underlying partonic cross sectionsσ. At hadron colliders, the LO production of W ± and Z bosons involves the annihilation of a quark-antiquark pair of the same (qq → Z + X) or different (qq → W + X) flavour. At NLO, the Born terms are supplemented with initial-state real gluon emission, virtual gluon exchange, and contributions from gluon-quark and gluonantiquark scattering processes. At NNLO, additional gluon radiation and virtual exchanges further contribute to the total cross section [15,18]. Although at LO the partonic cross sections are independent of α S , the vertices of the higher-order terms introduce a dependence on α S that enables the determination of the QCD coupling by comparing high-precision theoretical calculations to the experimental data. The size of such higher-order corrections [19], encoded in the so-called K-factor, amounts to K = σ nnlo /σ lo ≈ 1.25-1.37 as derived with mcfm v.8.0 [16] for the W ± and Z cross sections measured at 7 and 8 TeV in the CMS fiducial acceptance, and indicates that EW boson production is indeed sensitive to α S (m Z ) at NNLO accuracy. In this work, the NNLO cross sections are computed with the mcfm code interfaced with lhapdf v.6.1.6 [22] to access four different PDFs: CT14 [23], HERAPDF2.0 [24], MMHT14 [25], and NNPDF3.0 [26]. All these PDFs use as the default central set the one with the QCD coupling constant fixed to α S (m Z ) = 0.118 in their global fits of the data, but also provide a variety of alternative sets with their corresponding central values derived for different fixed values of α S (m Z ). We note, however, that when the QCD coupling constant is left free in their NNLO PDF fits, the following values are preferred by the different PDF sets: α S (m Z ) = 0.115 (CT14) [23], 0.108 (HERAPDF2.0) [24], and 0.1172 (MMHT2014) [25]. The HERAPDF2.0 set is obtained from fits to HERA deep inelastic scattering (DIS) data only. The CT14, MMHT14, and NNPDF3.0 global fits have been obtained including DIS, fixed target, and LHC measurements. These latter PDF sets incorporate one or two W ± or Z differential CMS distributions at 7 TeV [13] in their global fits, but did not use any of the twelve absolute inclusive EW boson cross sections listed in table 1, and therefore the corresponding values of α S extracted here are truly independent of the data contributing to the extraction of PDF sets themselves. The so-called G µ electroweak scheme, where the input parameters are m W , m Z , and G F , is used in all the predictions. The leptonic W and Z branching fractions are obtained in mcfm from the theoretical leptonic width (computed at LO in electroweak accuracy) normalized to the total W and Z widths experimentally measured [1]. All numerical results have been obtained using the latest SM parameters for particle masses, widths, and couplings [1]. For simplicity, the default value of the charm quark mass in mcfm and NNPDF3.0, m c = 1.275 GeV,

JHEP06(2020)018
is used for all the calculations -rather than the preferred values of the other PDF sets: m c = 1.3 GeV (CT14), 1.43 GeV (HERAPDF2.0), 1.4 GeV (MMHT14) -because the computed cross sections vary only by a few per mille, within the mcfm numerical uncertainty.
The default renormalisation and factorisation scales are set to the corresponding EW boson mass for each process, µ = µ f = µ r = m W , m Z . For all PDF sets, we computed the NNLO cross sections at various α S (m Z ) values over the range [0.115-0.121]. For NNPDF3.0, the available values are α S (m Z ) = 0. 115, 0.117, 0.118, 0.119, and 0.121; and for the other PDF sets they are α S (m Z ) = 0. 115, 0.116, 0.117, 0.118, 0.119, 0.120, and 0.121. Technically, the central sets selected via lhapdf for this study are: CT14nnlo as 0iii (for iii = 115-121), HERAPDF20 NNLO ALPHAS iii (for iii = 115-121), MMHT2014nnlo asmzlargerange (with α S (m Z ) = 0.115, . . . , 0.121 grids), and NNPDF30 nnlo as 0iii (with iii = 115-121). For the PDF uncertainties, only NNPDF3.0 provides independent replicas for each α S (m Z ) set, which we use in our calculations and uncertainties propagation, whereas the rest of PDFs use the same eigenvalues corresponding to the set determined with α S (m Z ) = 0.118. Calculations are carried out implementing the fiducial selection criteria for the final-state charged leptons corresponding to each of the six different measurements (W + e , W − e , Z e , W + µ , W − µ , Z µ ) at √ s = 7 and 8 TeV listed in table 1, thereby providing altogether twelve theoretical cross sections per PDF that can be used to individually extract α S (m Z ).
The PDF uncertainties of the theoretical fiducial cross sections are obtained by taking into account the different eigenvector sets, or replicas, that come with each of the PDFs. We use the "official" prescriptions of each PDF set to compute their associated uncertainties. More specifically, the PDF uncertainties are calculated from the cross sections obtained with the central PDF member (σ 0 ) and with the rest of eigenvalues or replicas (σ i ) as follows: • For CT14, the uncertainty eigenvectors are considered in pairs from the i = 1-56 PDF members. The largest positive and negative differences from each pair are summed quadratically to obtain the corresponding positive and negative PDF uncertainties: The CT14 PDF set results in asymmetric uncertainties interpreted as a 90% confidence level interval. To convert those to one standard deviation, as for the rest of PDF sets, they are divided by a factor of √ 2 erf −1 (0.9) ≈ 1.645.
• For HERAPDF2.0, a first asymmetric uncertainty is derived from the so-called 'EIG' (experimental uncertainties) PDF members, and a second one from the i = 1-10 'VAR' (variation) members, as for CT14. A third asymmetric uncertainty is taken from the i = 11-13 VAR members, as the maximum positive and negative differences ∆σ i with respect to σ 0 . Finally, all positive and negative uncertainties are separately added quadratically to get the final uncertainties.
To determine the scale uncertainty associated with missing corrections beyond the NNLO accuracy, the mcfm cross sections are recalculated for each PDF and measurement using factorisation and renormalisation scales varied within factors of two, such that the ratio of the two scales is not less than 0.5 or more than 2. This gives seven combinations: The largest and smallest cross sections of the seven combinations are determined, and the scale uncertainty is taken as half the difference of the extremal values. The scale variation uncertainties amount to 0.5-1% of the theoretical cross sections.
Since mcfm does not include EW corrections, arising from additional W ± , Z, and/or photons exchanged and/or radiated in the partonic process, those are computed separately. For this purpose, the mcsanc v.1.01 code [17] is used. For e ± final states, we follow the "calorimetric" prescription, proposed in theoretical W ± and Z boson production cross section benchmarking studies at the LHC [27], and recombine any radiated photon with the e ± if their relative distance in the pseudorapidity-azimuth plane is For µ ± final states, we use directly the "bare" mcsanc cross section. We run mcsanc at NLO with EW corrections on and off, and compute the corresponding multiplicative factor K EW = σ(nlo,ew on)/σ(nlo,ew off), which is used to correct the pQCD mcfm results. The EW corrections, in the range of 1-4%, are all negative, i.e. they reduce the overall cross section with respect to the pure pQCD result. Since the EW corrections are small, their associated uncertainties are neglected hereafter because they would propagate into the final computed W ± and Z cross section at a few per mille level, below the numerical uncertainty of the mcfm calculation. Subtracting the EW corrections, rather than applying them multiplicatively via a K EW factor as done here, gives consistent results within the (neglected) per mille uncertainties. For simplicity, in the mcsanc calculations, the electron pseudorapidity range 1.44 < |η| < 1.57 (excluded in the actual measurements) is also included, since we are interested in the relative correction, this small range (present in both the numerator and denominator of the correcting factor) does not affect the K EW ratio. The roles of photon-induced contributions and of mixed QCD⊕QED NLO corrections to Drell-Yan processes in pp collisions have been computed in refs. [28] and [29], respectively. The impact of such corrections to the inclusive W ± and Z cross sections is at a few per mille level, and also neglected here.
All the relevant sources of uncertainties in the W ± and Z boson cross sections are summarised in table 2. The largest experimental and theoretical uncertainties come from the integrated luminosity and PDF knowledge, respectively. Each calculated cross section has a numerical accuracy, as reported by mcfm, in the range of 0.2-0.6%. Such an uncertainty is commensurate with the typical differences encountered when computing NNLO W ± and -6 -JHEP06(2020)018 Z boson cross sections with different pQCD codes that implement higher-order virtual-real corrections with various methods [27].

EW boson fiducial cross sections: data versus theory
All the experimental and theoretical fiducial cross sections for W + , W − , and Z production in pp collisions are given in tables 3 and 4 for 7 and 8 TeV, respectively. For each measurement, the fiducial cross section definition and the experimental result are listed along with their uncertainties from the different sources listed in table 2. The theoretical mcfm predictions computed with all four PDF sets for their preferred default α S (m Z ) = 0.118 value are listed including their associated PDF, α S (obtained, as described in section 5.2, from the cross section change when α S (m Z ) is modified by ±0.001), and scale uncertainties.
For each system, the NLO mcsanc EW correction factors (absolute and relative) are also listed. For the results at √ s = 8 TeV, the theoretical result obtained with the alternative fewz NNLO pQCD calculator [30], using the MSTW2008 PDF set [31] as provided in the original ref. [14], is also listed to show the very similar theoretical predictions expected with an alternative NNLO code and a pre-LHC PDF set. For each of the twelve experimental W ± and Z boson cross section measurements listed in table 1, we have computed their corresponding theoretical NNLO pQCD predictions using the four PDF sets and five to seven different values of α S (m Z ). It is important to stress again that, for each QCD coupling constant, we use the specific PDF sets that are associated with that particular α S (m Z ) value. We calculated the NLO EW corrections using NNPDF3.0 and α S (m Z ) = 0.118. By comparing the whole set of theoretical calculations to the experimental data, a preferred value of the QCD coupling constant can be derived for each PDF set as explained in the next section.

System
Fiducial cross section Table 3. Experimental and theoretical fiducial cross sections for W ± and Z production in pp collisions at √ s = 7 TeV, with the uncertainty sources listed in table 2. The NNLO pQCD results are obtained with mcfm for α S (m Z ) = 0.118 using the CT14, HERAPDF2.0, MMHT14, and NNPDF3.0 PDF sets. (The quoted α S uncertainties are derived from the cross section changes when α S (m Z ) is modified by ±0.001). The NLO EW corrections are computed with mcsanc.
ellipses represent the contours of the joint probability density functions (Jpdfs) of the theoretical and experimental results, with a width representing a two-dimensional one standard deviation obtained from the product of the probability densities of the experimental and -8 -

System
Fiducial cross section 2240 ± 10 (stat) ± 20 (syst) ± 60 (lumi) pb 400 ± 10 (stat) ± 10 (syst) ± 10 (lumi) pb Table 4. Experimental and theoretical fiducial cross sections for W ± and Z production in pp collisions at √ s = 8 TeV, with the uncertainty sources listed in table 2. The NNLO pQCD results are obtained with mcfm for α S (m Z ) = 0.118 using the CT14, HERAPDF2.0, MMHT14, and NNPDF3.0 PDF sets, as well as with fewz using the MSTW2008 PDF. (The quoted α S uncertainties are derived from the cross section changes when α S (m Z ) is modified by ±0.001). The NLO EW corrections are computed with mcsanc. predictions tending to be systematically above the data, and the NNPDF3.0 ones below the latter. In between the results of these two PDF sets, the cross sections derived with MMHT14 tend to be above those with CT14, although they are often very similar and overlap most of the time. Alternatively, the results of figures 1-2 indicate that, in order to reproduce the experimental cross sections, HERAPDF2.0 (NNPDF3.0) tends in general to prefer a smaller (larger) value of α S (m Z ) than other PDFs, and that the predictions from CT14 and MMHT14 tend to be less scattered over the α S (m Z ) axis than those from HERAPDF2.0 and NNPDF3.0. The HERAPDF2.0 (MMHT14) filled ellipses have the smallest (largest) relative slope as a function of α S (m Z  Table 5. Overall goodness-of-fit per number of degrees of freedom, χ 2 /ndf, among the twelve experimental measurements of W ± and Z boson production cross sections and the corresponding theoretical calculations obtained with the four different PDF sets for their default α S (m Z ) = 0.118 value. The first (second) row is obtained symmetrising the PDF uncertainties of the cross sections obtained with the CT14, HERAPDF2.0, and MMHT14 sets to the largest (smallest) of their respective values.
extracting the strong coupling constant, because this means that the underlying α S (m Z ) value in the calculations has a larger impact on the computed cross sections, also leading to a lower propagated uncertainty in the α S (m Z ) value derived by comparing the theoretical prediction to the experimental value. Overall, the theoretical predictions computed using the world-average value of the QCD coupling constant (vertical dashed line in figures 1-2) agree well with the experimental values within the uncertainties. The level of data-theory agreement can be quantified with a goodness-of-fit test, χ 2 = ξ i (M −1 ) ij ξ j , where M is the covariance matrix taking all the uncertainties and their correlations into account, as explained in section 5.3, and ξ i = σ i,th − σ i,exp is the difference between theoretical and experimental cross sections for each PDF set. In the χ 2 calculation, the asymmetric uncertainties of the CT14, HERAPDF2.0, and MMHT14 PDF sets are symmetrised to the largest of the two values and also separately to the smallest of the two values. The results are listed in table 5.
5 Extraction of the QCD coupling constant 5.1 Extraction of α S (m Z ) for each single W ± and Z cross section measurement The dependence of the theoretical cross sections on the QCD coupling constant α S (m Z ), shown in figures 1-2, is fitted through a linear χ 2 -minimisation procedure over α S (m Z ) ∈ [0.115, 0.121], to extract the slope k. Over the considered α S (m Z ) range, the empirical linear fit describes well the observed α S -dependence of the theoretical cross section for all PDF sets. The value of α S (m Z ) preferred by each individual measurement is determined by the crossing point of the fitted linear theoretical curve with the experimental horizontal line. The resulting α S (m Z ) values are listed in table 6. For each theoretical point used in the fit, the uncertainty in the cross section is given by the quadratic sum of its associated PDF, scale, and numerical uncertainties. In order to exploit the dependence of σ th (V) on α S to quantitatively derive the latter, a joint probability density function is constructed for each PDF prediction, as explained next. When the theoretical cross section value has positive and negative uncertainties δ +-, the experimental cross section is σ exp with uncertainty δ exp , the slope of the fit is k, and the fitted QCD coupling constant is α S (m Z ), then the Jpdf as a function of σ and α S is proportional to  Table 6. Extracted α S (m Z ) values from the different data-theory W ± and Z boson production cross section comparisons for each PDF set, with associated uncertainties from different experimental (statistical, integrated luminosity, and systematic) and theoretical (PDF, scale, and numerical) sources.

JHEP06(2020)018
The sign of (σ − σ exp − k(α S − α S (m Z ))) determines which of the δ +-is used. For symmetric uncertainties the Jpdfs have elliptical contours, but for asymmetric ones they are two filled ellipses combined together. This procedure is repeated for all the twelve different measurements and for all four PDF sets, and plotted as the filled ellipses shown in figures 1-2, where each coloured area corresponds to one two-dimensional (in σ and α S (m Z )) standard deviation.

Propagation of α S (m Z ) uncertainties
Appropriate propagation of the separated experimental and theoretical uncertainties into each value of α S (m Z ) obtained from each particular W ± and Z measurement, is crucial to combine all estimates taking into account their correlations, and extract a single final α S (m Z ) result. The method employed to determine the individual sources of uncertainties associated with a given α S (m Z ) value is similar to that used in refs. [8,9] for the α S (m Z ) extraction from inclusive tt cross sections. In summary, each source of uncertainty δσ propagates into a corresponding α S (m Z ) uncertainty through δσ/k, where k is the slope of the fit of the theoretical cross section versus α S (m Z ). The validity of such a simple propagation of uncertainties can be demonstrated calculating first a theoretical uncertainty distribution by adding up quadratically the PDF, scale, and numerical uncertainties assuming Gaussian distributions (for the asymmetric PDF uncertainties, only the largest of the positive and negative uncertainties are used). Then, a Jpdf can be derived from the product of the theoretical and experimental distributions f th (σ|α S (m Z )) and f exp (σ). Integration over σ gives the marginalised posterior The expected value of the theoretical probability distribution changes linearly according to the fitted first-order polynomial, but all the theoretical and experimental uncertainties remain the same for all α S (m Z ) values. Since all the uncertainties have Gaussian distributions and the marginalisation is, in essence, a convolution, the resulting α S (m Z ) posterior will be also Gaussian, with the impact of each cross section uncertainty adding quadratically to the α S (m Z ) uncertainty. More specifically, each cross section uncertainty source δσ will result in a propagated α S (m Z ) uncertainty in δσ/k size, where k is the slope of the linear fit to theoretical calculations. In this demonstration, we symmetrised the PDF uncertainties for simplicity, but the δσ/k prescription will be also used hereafter for the case of asymmetric uncertainties. All the extracted α S (m Z ) values, along with the uncertainty breakdowns from every source, for each system and PDF set are given in table 6. The results from the MMHT14 PDF feature the extracted α S (m Z ) values with the lowest overall uncertainty, in some cases as low as 3%.

Combination of all individual α S (m Z ) values
From the twelve α S (m Z ) extractions per PDF set listed in table 6, we can determine a single α S (m Z ) value by appropriately combining them taking into account their uncorrelated, partially-, and fully-correlated experimental and theoretical uncertainties. For this, -14 -

JHEP06(2020)018
the program convino v1.2 [32] is employed, which uses a χ 2 minimisation to determine the best estimate. In the current analysis, the Neyman χ 2 code option is always used. As an independent cross-check, we confirm that, for symmetric uncertainties, the results are identical to those obtained with the BLUE method [33]. The following correlation coefficients are used: • The integrated luminosity uncertainty is fully correlated for all α S (m Z ) results obtained at the same √ s, but fully uncorrelated between the two different c.m. energies.
• The experimental systematic uncertainty is partially correlated. Since the exact correlation values impact the final result, a dedicated study of the correlations is carried out in section 5.3.1.
• The experimental statistical uncertainty is fully uncorrelated among α S (m Z ) extractions.
• The PDF uncertainty is partially correlated for the α S (m Z ) values extracted with the same PDF set, as discussed in detail in section 5.3.2.
• The scale uncertainty is partially correlated, as explained in section 5.3.2.
• The theoretical numerical uncertainty is fully uncorrelated among α S (m Z ) extractions.
By properly implementing all the uncertainties and their correlations in convino, we can derive a single final combined α S (m Z ) value and associated uncertainties for each PDF set.

Correlations among the experimental systematic uncertainties
For all the experimental measurements of the cross sections, the size of their systematic uncertainties of each type are listed in table 7. The absolute uncertainties are given in the same proportions as in table 7, but rescaled such that they add up quadratically to the total experimental systematic uncertainty listed in tables 3 and 4. The detailed correlations between the different uncertainty sources of table 7 are listed in tables 10 to 14 of the appendix. By using the experimental systematic uncertainty breakdown and the correlations between the uncertainties, the total correlations between the experimental uncertainty sources can be calculated using the formula where the subscript k labels the uncertainty (e.g. background subtraction), and i, j denote the associated measurement (e.g. W + e at 7 TeV). The calculated total correlations among experimental systematic uncertainties are given in table 15 of the appendix. Many of the propagated experimental uncertainties appear strongly correlated. To give an idea of the correlations among α S (m Z ) estimates taking into account all the uncertainty sources, both experimental and theoretical, an example correlation matrix for the NNPDF3.0 set is given in table 16 of the appendix. One can see that many of the α S (m Z ) values derived within a given PDF set are strongly correlated, especially across the same √ s.

Correlations among PDF and scale uncertainties
In the theoretical cross section calculations, the PDF uncertainties are in the range of a few percent, scale uncertainty up to one percent, and numerical uncertainty around half a percent (figures 3 and 4). The mcfm numerical uncertainty cannot be neglected because differences with respect to the prediction computed with the central eigenvalue/replica members are used to calculate the PDF uncertainties. The cross sections for the central PDF members have intrinsic numerical fluctuations that impact the PDF uncertainty magnitude, asymmetry, and also the correlations among theoretical uncertainties. To take this into account, a Pearson correlation coefficient is calculated for each pair of measurements and for each PDF set using the cross sections from all the PDF members that were used in calculating the PDF uncertainty. Similarly for the scale uncertainties, for each pair of measurements the Pearson correlation coefficient is calculated using the results obtained from varying the theoretical scales. The correlations are mostly in the 0. 8-0.9, 0.4-0.7, 0.2-0.6, and 0.9-1.0 ranges for CT14, HERAPDF2.0, MMHT14, and NNPDF3.0, respectively. The scale correlations are around 0.6-0.9. When combining the α S (m Z ) estimates, the specific correlation coefficient calculated for every specific pair of estimates is used.     Table 9. Sensitivity of the final α S (m Z ) extractions per PDF set to various data, uncertainties, and correlation assumptions. Top rows: extractions of α S (m Z ) using only the 7 and 8 TeV measurements separately. Bottom rows: extractions of α S (m Z ) when symmetrising the asymmetric PDF uncertainties by taking the maximum of the negative and positive values (left), and when adding a 1% uncorrelated uncertainty to all cross sections (right).

Results and discussion
coupling constant derived with HERAPDF2.0, α S (m Z ) = 0.1072 +0.0043 −0.0040 , is between 1.7 and 2.7 standard deviations smaller than the rest of extractions (figure 4 left). Although as discussed later, in the context of the cross-checks described in table 9, such a disagreement is reduced when symmetrising the HERAPDF2.0 uncertainties to their maximum values. As discussed before, since the HERAPDF2.0 cross sections for α S (m Z ) = 0.118 tend to overpredict the measured W ± and Z boson production cross sections (figures 1-2), a datatheory agreement can only be obtained for a value of α S that is reduced compared to its -18 -JHEP06(2020)018 default value. For all global PDF fits extracted with different α S values as input, there exists a generic anticorrelation between the values of α S (Q 2 ) and the parton densities evaluated at (x, Q 2 ), particularly, for the gluon and in turn (through perturbative evolution) for the sea quarks. It is thereby important to analyze in more detail the differences between the PDF sets for each flavour. For this purpose, a comparison study of parton densities and parton luminosities has been carried out with apfel v2.7.1 [34]. This study indicates that the HERAPDF2.0 u-quark densities (and the overall quark-antiquark luminosities) are enhanced by about 5% compared to the rest of PDFs in the (x, Q 2 ) region of relevance for EW boson production. This fact increases the weight of the LO contributions to the theoretical W ± and Z boson production cross sections, and thereby pushes down the cross section contributions from higher-order pQCD diagrams that are sensitive to α S (m Z ). The effective result is a comparatively reduced α S (m Z ) value. The level of agreement between the twelve individual and the total α S (m Z ) extractions turns out to be good for this PDF set (χ 2 /ndf ≈ 1 in table 8), because of the relatively wide span of α S values derived and their associated large uncertainties (figure 3). Since HERAPDF2.0 uses DIS data alone, and therefore lacks the extra constraints on the PDFs provided by the LHC data, we conclude that one would need an updated refit of these parton densities to an extended set of experimental data, including LHC results, before relying on the QCD coupling constant derived following the procedure described here.
The features of the α S (m Z ) results obtained with NNPDF3.0 show the opposite behaviour to those observed for the HERAPDF2.0 set. The W ± and Z boson production cross sections computed with this PDF tend to underpredict the experimental measurements (figures 1-2), and yield an overall bad data-theory agreement (χ 2 /ndf ≈ 3 in table 5), for the baseline α S (m Z ) = 0.118 coupling constant. A reproduction of the individual measurements by theory can thus be achieved only for an α S (m Z ) value that is enlarged compared to the default value for this PDF set. Thus, many of the individual α S (m Z ) extractions obtained with NNPDF3.0 have values relatively larger than those obtained for other PDFs. However, the final combined NNPDF3.0 value appears shifted down to α S (m Z ) = 0.1147 ± 0.0023 (table 8), falling outside of the region around α S (m Z ) ≈ 0.120 defined by most of the individual estimates (table 6), and, consequently, the final level of agreement of the combined and single extractions is poor (χ 2 /ndf ≈ 3 in table 8). Such a seemingly counterintuitive behaviour is due to the presence of strong correlations among individual extractions for a fixed value of α S (m Z ), and the fact that the lowest α S (m Z ) values derived have smaller uncertainties than the rest, thereby pulling down the final average. The effect of a combined result lying outside of the range of the input values caused by, e.g. large underlying nonlinearities among individual estimates, is called "Peelle's pertinent puzzle" [35]. The absence of parametrisation bias in this neural-network PDF results in parton densities that can be less well constrained (have larger uncertainties) than the rest of PDFs in some regions of phase space. The same apfel v2.7.1 study mentioned above indicates that the NNPDF3.0 quark-antiquark luminosities tend to be somewhat less precise than those from the other PDF sets in the relevant (x, Q 2 ) range for W ± and Z production. A larger span of replicas results in nontrivial correlations among PDF uncertainties that push the final α S (m Z ) value off the individual extractions for each single -19 -JHEP06(2020)018 measurement. In any case, the latest version released (v3.1) of the NNPDF global fit [36] shows much better agreement of the theoretical EW boson cross sections with the LHC data for the central α S (m Z ) = 0.118 value. However, this latter NNPDF3.1 set cannot be directly employed to independently extract α S (m Z ) from the CMS measurements through the approach discussed here, as this updated global fit does include already the absolute normalisation of a fraction of the W ± and Z cross sections used in this work.
The final α S (m Z ) extractions are plotted in figure 4 (left) -with (asymmetric, where needed) parabolas constructed so as to have a minimum at each final central α S (m Z ) result and (one standard deviation) uncertainties matching those listed in table 8 -compared with the current world average of α S (m Z ) = 0.1181 ± 0.0011 (orange band).
To analyse the robustness and stability of the final α S (m Z ) extractions to the underlying data sets, their uncertainties, and correlations, we repeat the convino combination varying several ingredients, as explained next. First, α S (m Z ) values are extracted using separately the measurements at √ s = 7 and 8 TeV alone, as shown in the top rows of table 9. This separation of data sets yields final α S (m Z ) values mostly consistent with those derived from the combined ones listed in table 8, with the largest deviations from the original results being those from the 7 TeV NNPDF3.0 and 8 TeV HERAPDF2.0 extractions, with 1.0 and 0.85 standard deviations, respectively. Although the CMS luminosity studies confirm that these uncertainties are fully uncorrelated between 7 and 8 TeV, we have checked the impact of relaxing such an assumption by assuming a 0.5 correlation factor between them. Such a correlation factor results in a change of the final α S (m Z ) by at most one-third of the current total uncertainty. Another cross-check is carried out by symmetrising the PDF uncertainties to their maximum value of the two (this does not apply to NNPDF3.0, because its uncertainties are symmetric by construction). The corresponding results are given in the left bottom half of table 9. Changing from asymmetric to symmetric PDF uncertainties causes the HERAPDF2.0 combined value to increase by 1.1 standard deviation, whereas all other PDF results are consistent with the default α S (m Z ) extractions. Such a large sensitivity to changes in the PDF uncertainty confirms the relative lack of robustness of the α S (m Z ) values derived for HERAPDF2.0 in our analysis, because the asymmetries of the PDF uncertainties can be significantly affected by random numerical errors. To further test the sensitivity of the α S (m Z ) extraction to the assumptions made on the underlying W ± and Z cross section uncertainties and their correlations, the original analysis is repeated by adding an uncorrelated 1% numerical uncertainty to all theoretical cross sections. Such a value accounts for possible overlooked small uncorrelated uncertainties, e.g. coming from the use of different codes for the theoretical pQCD and/or EW calculations [27]. The impact of such a change is not significant in the final results, as observed by comparing the numbers in table 8 and those in the bottom-right columns of table 9. Further similar tests and cross-checks have been carried out in a recent α S (m Z ) extraction that exploits all the LHC electroweak boson data [40] with the same approach used here, which confirm these conclusions. All these systematic tests indicate that HERAPDF2.0 and NNPDF3.0 have larger variations when changing the ingredients of the combination, but for CT14 and MMHT14 the final α S (m Z ) values extracted are reasonably robust within the quoted uncertainties.

JHEP06(2020)018
Among PDFs, the results obtained using MMHT14 and CT14 feature the largest sensitivity to α S variations, i.e. they show a larger k slope, eq. (5.1), compared to those obtained with HERAPDF2.0 and NNPDF3.0 (figures [1][2]. Since the uncertainty in the α S (m Z ) value derived from HERAPDF2.0 is the largest (up to twice as large as some of the other extractions), because of the absence of constraining LHC input data in this HERA-only PDF fit, and since the final NNPDF3.0 result has a larger tension between the combined and individual extractions from each single measurement (table 8), we consider the values extracted with CT14, α S (m Z ) = 0.1163 +0.0024 −0.0031 , and MMHT14, α S (m Z ) = 0.1186 ± 0.0025, as the most reliable in this analysis. Providing a single final α S (m Z ) value from this study is not obvious because, in general, there is no unique way to derive a final best estimate of α S based on the results obtained from different PDF sets. An unbiased approach for combining results from different PDFs, in line with the PDG practice [1] as well as with the procedure employed to produce the PDF4LHC combined PDF set [37], is to average them without applying any further weighting. The same approach was followed also in the similar combination of QCD coupling constant values obtained from the inclusive tt cross sections [9]. By taking the straight average of the mean values and of the uncertainties of the individual CT14 and MMHT14 combinations, we obtain a final value of the QCD coupling constant at the Z pole mass, α S (m Z ) = 0.1175 +0.0025 −0.0028 , with a total (symmetrised) uncertainty of 2.3%. Such a result compares very well with the α S (m Z ) = 0.1177 +0.0034 −0.0036 value, with an uncertainty of ≈3%, extracted from the theoretical analysis of top pair cross section data [9]. The right plot of figure 4 shows the α S (m Z ) parabola extracted combining the CT14 and MMHT14 results. This final extraction is fully consistent with the PDG world average (orange band), and has an overall uncertainty similar to that of other recent determinations at this level of (NNLO) theoretical accuracy, such as those from EW precision fits [38], and tt cross sections [8][9][10].

Summary
We have used twelve measurements of the inclusive fiducial W ± and Z production cross sections in proton-proton collisions (pp) at √ s = 7 and 8 TeV, carried out in the electron and muon decay channels by the CMS experiment, to extract the value of the strong coupling constant at the Z pole mass, α S (m Z ). The procedure is based on a detailed comparison of the measured electroweak boson cross sections to theoretical calculations computed at next-to-next-to-leading-order accuracy with the CT14, HERAPDF2.0, MMHT14, and NNPDF3.0 parton distribution function (PDF) sets. The overall data-theory agreement is good within the experimental and theoretical uncertainties. A χ 2 -minimisation procedure has been employed to combine all twelve individual α S extractions per PDF set, properly taking into account all individual sources of experimental and theoretical uncertainties, and their correlations. The following combined values are extracted for the four different PDFs: α S (m Z ) = 0.1163 +0.0024 −0.0031 (CT14), 0.1072 +0.0043 −0.0040 (HERAPDF2.0), 0.1186 ± 0.0025 (MMHT14), and 0.1147 ± 0.0023 (NNPDF3.0). The largest propagated uncertainties are associated with the experimental integrated luminosity and theoretical intra-PDF uncertainties. Among the four extractions, the cross section calculated with the CT14 and -21 -

JHEP06(2020)018
MMHT14 sets appear as the most sensitive to the underlying α S value and, at the same time, the derived α S (m Z ) values are the most robust and stable with respect to variations in the data and theoretical cross sections, their uncertainties, and correlations. The result derived combining the CT14 and MMHT14 extractions, α S (m Z ) = 0.1175 +0.0025 −0.0028 , has a 2.3% uncertainty that is comparable to that previously obtained in a similar analysis of the inclusive tt cross sections in pp collisions at the LHC. This extracted value is fully compatible with the current α S (m Z ) world average.

Acknowledgments
We congratulate our colleagues in the CERN accelerator departments for the excellent performance of the LHC and thank the technical and administrative staffs at CERN and at other CMS institutes for their contributions to the success of the CMS effort. In addition, we gratefully acknowledge the computing centres and personnel of the Worldwide LHC Computing Grid for delivering so effectively the computing infrastructure essential to our analyses. Finally, we acknowledge the enduring support for the construction and operation of the LHC and the CMS detector provided by the following funding agencies: BMBWF and FWF (Austria

A Correlation matrices of the experimental measurements
Relevant correlation matrices of the experimental systematic uncertainties in the measurements of W ± and Z boson production cross sections, and in their associated extractions of α S (m Z ), are listed in tables 10-14 and 15-16, respectively.            Table 16. Total experimental and theoretical correlations between all the α S (m Z ) values extracted from all measurements at 7 and 8 TeV [13,14] of W ± and Z boson production cross sections, obtained using the NNPDF3.0 PDF set.

JHEP06(2020)018
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.
[13] CMS collaboration, Measurement of the inclusive W and Z production cross sections in pp collisions at √ s = 7 TeV, JHEP 10 ( -27 -