Measurement of $\mathrm{t\bar t}$ normalised multi-differential cross sections in pp collisions at $\sqrt s=13$ TeV, and simultaneous determination of the strong coupling strength, top quark pole mass, and parton distribution functions

Normalised multi-differential cross sections for top quark pair ($\mathrm{t\overline{t}}$) production are measured in proton-proton collisions at a centre-of-mass energy of 13 TeV using events containing two oppositely charged leptons. The analysed data were recorded with the CMS detector in 2016 and correspond to an integrated luminosity of 35.9 fb$^{-1}$. The double-differential $\mathrm{t\overline{t}}$ cross section is measured as a function of the kinematic properties of the top quark and of the $\mathrm{t\overline{t}}$ system at parton level in the full phase space. A triple-differential measurement is performed as a function of the invariant mass and rapidity of the $\mathrm{t\overline{t}}$ system and the multiplicity of additional jets at particle level. The data are compared to predictions of Monte Carlo event generators that complement next-to-leading-order (NLO) quantum chromodynamics (QCD) calculations with parton showers. Together with a fixed-order NLO QCD calculation, the triple-differential measurement is used to extract values of the strong coupling strength $\alpha_S$ and the top quark pole mass ($m_\mathrm{T}^\text{pole}$) using several sets of parton distribution functions (PDFs). Furthermore, a simultaneous fit of the PDFs, $\alpha_S$, and $m_\mathrm{T}^\text{pole}$ is performed at NLO, demonstrating that the new data have significant impact on the gluon PDF, and at the same time allow an accurate determination of $\alpha_S$ and $m_\mathrm{T}^\text{pole}$.


Introduction
Measurements of top quark pair (tt) production are important for checking the validity of the standard model (SM) and searching for new phenomena. In particular, the large data set delivered by the CERN LHC allows precise measurements of the tt production cross section as a function of tt kinematic observables. These can be used to check the most recent predictions of perturbative quantum chromodynamics (QCD) and to constrain input parameters, some of which are fundamental to the SM. At the LHC, top quarks are predominantly produced via gluon-gluon fusion. Using measurements of the production cross section in a global fit of parton distribution functions (PDFs) can help determine the gluon distribution at large values of x [1][2][3], where x is the fraction of the proton momentum carried by a parton. Furthermore, measurements of the cross section as a function of the tt invariant mass, from the threshold to the TeV region, provide high sensitivity for constraining the top quark pole mass, m pole t , which is defined as the pole of the top quark propagator (see e.g. Refs. [4][5][6]). At LHC energies, a large fraction of tt events is produced with additional hard jets in the final state. Events containing such additional jets constitute important backgrounds for interesting but rare SM processes such as the associated production of a Higgs boson and tt, as well as for searches for new physics associated with tt production, and must therefore be well understood. Within the SM, processes with extra jets can also be used to constrain the strong coupling strength, α S , at the scale of the top quark mass. Furthermore, the production of tt in association with extra jets provides additional sensitivity to m pole t since gluon radiation depends on m pole t through threshold and cone effects [7]. Differential cross sections for tt production have been measured previously in protonantiproton collisions at the Tevatron at a centre-of-mass energy of 1.96 TeV [8,9] and in protonproton (pp) collisions at the LHC at √ s = 7 TeV [10-14], [8][9][10][11][12][13][14][15][16][17][18][19][20], and 13 TeV [21][22][23][24][25][26]. A milestone was reached in three CMS analyses [20][21][22], where the tt production dynamics was probed with double-differential cross sections. The first analysis [20] used data recorded at √ s = 8 TeV by the CMS experiment in 2012. Only tt decays where, after the decay of each top quark into a bottom quark and a W boson, both of the W bosons decay leptonically were considered. Specifically, the e ± µ ∓ decay mode (eµ) was selected, requiring two oppositely charged leptons and at least two jets. Our present paper provides a new measurement, following the procedures of Ref. [20]. It is based on data taken by the CMS experiment in 2016 at √ s = 13 TeV, corresponding to an integrated luminosity of 35.9 ± 0.9 fb −1 . In addition to eµ, the decay modes e + e − (ee) and µ + µ − (µµ) are also selected, roughly doubling thereby the total number of expected tt signal events. Our latest measurement complements the analyses [21,22], based on data taken at √ s = 13 TeV, but using tt decays in the +jets ( = e, µ) final state.
As in the previous work [20], measurements are performed of the normalised double-differential tt production cross section as a function of observables describing the kinematic properties of the top quark, top antiquark, and the tt system: the transverse momentum of the top quark, p T (t), rapidity of the top quark, y(t); the transverse momentum, p T (tt), the rapidity, y(tt), and invariant mass, M(tt), of the tt system; the pseudorapidity difference between the top quark and antiquark, ∆η(t, t), and the angle between the top quark and antiquark in the transverse plane, ∆φ(t, t). When referring to the kinematic variables p T (t) and y(t), we use only the parameters of the top quark and not of the top antiquark, to avoid double counting of events. In all, the double-differential tt cross section is measured as a function of six different pairs of kinematic variables. As demonstrated in Ref. [20], the different combinations of kinematic variables are sensitive to different aspects of the QCD calculations.
For the first time at the LHC, the triple-differential cross section is measured as a function of M(tt), y(tt), and N jet , where N jet is the number of extra jets not arising from the decay of the tt system. For this purpose, the kinematic reconstruction algorithm is optimised to determine the invariant mass of the tt system in an unbiased way. As will be shown below, the tripledifferential measurements provide tight constraints on the parametrised gluon PDF, as well as on α S , and m pole t . Previous studies of additional jet activity in tt events at the LHC can be found in Refs. [22,27,28]. The α S and m pole t parameters were also extracted from measurements of the total inclusive tt production cross sections in Refs. [29][30][31][32][33][34].
The measurements are defined at parton level and must therefore be corrected for effects of hadronisation, and detector resolution and inefficiency. A regularised unfolding process is performed simultaneously in bins of the two or three variables in which the cross sections are measured. The normalised differential tt cross section is determined by dividing the distributions by the measured total inclusive tt production cross section, where the latter is evaluated by integrating over all bins in the respective observables.
The parton-level results are compared with theoretical predictions obtained with the generators POWHEG (version 2) [35,36] and MG5 aMC@NLO [37], interfaced to PYTHIA [38,39] for parton showering, hadronisation, and multiple-parton interactions (MPIs). They are also compared to theoretical predictions obtained at next-to-leading-order (NLO) QCD using several sets of PDFs, after applying corrections for non-perturbative (NP) effects.
The structure of the paper is as follows: Section 2 contains a brief description of the CMS detector. Details of the event simulation are given in Section 3. The event selection, kinematic reconstruction, and comparison between data and simulation are described in Section 4. The unfolding procedure is detailed in Section 5, the method to determine the differential cross sections is presented in Section 6, and the assessment of the systematic uncertainties is discussed in Section 7. We show the results of the measurement and their comparison to theoretical predictions in Section 8. Section 9 presents the extraction of α S and m pole t from the measured tt cross section, using several sets of PDFs, and Section 10 presents the simultaneous fit of the PDFs, α S , and m pole t to the data. Finally, Section 11 provides a summary.

The CMS detector
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 solenoid volume are a silicon pixel and strip tracker, a lead tungstate crystal electromagnetic calorimeter (ECAL), and a brass and scintillator hadron calorimeter, each composed of a barrel and two endcap sections. Forward calorimeters extend the η coverage provided by the barrel and endcap detectors. Muons are measured in gas-ionisation detectors embedded in the steel flux-return yoke outside the solenoid. Events of interest are selected using a two-tiered trigger system [40]. The first level, composed of custom hardware processors, uses information from the calorimeters and muon detectors to select events at a rate of around 100 kHz within a time interval of less than 4 µs. The second level, known as the high-level trigger (HLT), consists of a farm of processors running a version of the full event reconstruction software optimised for fast processing, and reduces the event rate to around 1 kHz before data storage. A more detailed description of the CMS detector, together with a definition of the coordinate system used and the relevant kinematic variables, can be found in Ref. [41].

Event simulation
Simulations of physics processes are performed with Monte Carlo (MC) event generators and serve three purposes: firstly, to obtain representative SM predictions of tt production cross sections to be compared to the results of this analysis. Secondly, when interfacing generated tt signal events with a detector simulation, to determine corrections for the effects of hadronisation, reconstruction and selection efficiencies, and resolutions that are to be applied to the data. Thirdly, when interfacing generated background processes to the detector simulation, to obtain predictions for the backgrounds. All MC programs used in this analysis perform the event generation in three subsequent steps: matrix-element (ME) level, parton showering matched to ME, and hadronisation. The tt signal processes are simulated with ME calculations at NLO in QCD. For all simulations the proton structure is described by the NNPDF 3.0 NLO PDF set with α S (m Z ) = 0.118 [42] where m Z = 91 GeV is the Zboson mass [43], and the value of the top quark mass is fixed to m MC t = 172.5 GeV. For the default signal simulation, the POWHEG (version 2) [35,44,45] generator is taken. The h damp parameter of POWHEG, which regulates the damping of real emissions in the NLO calculation when matching to the parton shower, is set to h damp = 1.581m MC t [46]. The PYTHIA program (version 8.2) [39] with the CUETP8M2T4 tune [46][47][48] is used to model parton showering, hadronisation and MPIs. An alternative sample is generated using the MG5 aMC@NLO (version 2.2.2) [37] generator, including up to two extra partons at the ME level at NLO precision. In this setup, referred to as MG5 aMC@NLO + PYTHIA, MADSPIN [49] is used to model the decays of the top quarks while preserving their spin correlation. The events are matched to PYTHIA using the FxFx prescription [50]. A second alternative sample is generated with POWHEG and interfaced with HERWIG++ (version 2.7.1) [51] using the EE5C tune [52].
The main background contributions originate from single top quarks produced in association with a W boson (tW), Z/γ * bosons produced with additional jets (Z+jets), W boson production with additional jets (W+jets) and diboson (WW, WZ, and ZZ) events. Other backgrounds are negligible. For all background samples, the NNPDF3.0 [42] PDF set is used and parton showering, hadronisation and MPIs are simulated with PYTHIA. Single top quark production is simulated with POWHEG (version 1) [36,53] using the CUETP8M2T4 tune in PYTHIA with the h damp parameter set to 172.5 GeV in POWHEG. The Z+jets process is simulated at NLO precision using MG5 aMC@NLO with up to two additional partons at ME level and matched to PYTHIA using the FxFx prescription. The W+jets process is simulated at leading order (LO) precision using MG5 aMC@NLO with up to four additional partons at ME level and matched to PYTHIA using the MLM prescription [54]. Diboson events are simulated with PYTHIA. Predictions are normalised based on their theoretical cross sections and the integrated luminosity of the data sample. The cross sections are calculated to approximate next-to-NLO (NNLO) for single top quark in the tW channel [55], NNLO for Z+jets and W+jets [56], and NLO for diboson production [57]. The tt simulation is normalised to a cross section of 832 +20 −29 (scale) ± 35(PDF+α S ) pb calculated with the TOP++ (version 2.0) program [58] at NNLO including resummation of nextto-next-to-leading-logarithm (NNLL) soft-gluon terms assuming m pole t = 172.5 GeV and the proton structure described by the CT14 NNLO PDF set [59].
To model the effect of additional pp interactions within the same bunch crossing (pileup), simulated minimum bias interactions are added to the simulated data. Events in the simulation are then weighted to reproduce the pileup distribution in the data, which is estimated from the measured bunch-to-bunch instantaneous luminosity assuming a total inelastic pp cross section of 69.2 mb [60].
In all cases, the interactions of particles with the CMS detector are simulated using GEANT4 (version 9.4) [61].

Event selection and tt kinematic reconstruction
The event selection procedure follows closely the one reported in Ref. [26]. Events are selected that correspond to the decay topology where both top quarks decay into a W boson and a b quark, and each of the W bosons decays directly into an electron or a muon and a neutrino. This defines the signal, while all other tt events, including those with at least one electron or muon originating from the decay of a τ lepton are regarded as background. The signal comprises three distinct final state channels: the same-flavour channels corresponding to two electrons (e + e − ) or two muons (µ + µ − ) and the different-flavour channel corresponding to one electron and one muon (e ± µ ∓ ). Final results are derived by combining the three channels.
At HLT level, events are selected either by single-lepton or dilepton triggers. The former require the presence of at least one electron or muon and the latter the presence of either two electrons, two muons, or an electron and a muon. For the single-electron and -muon triggers, p T thresholds of 27 and 24 GeV are applied, respectively. The same-flavour dilepton triggers require either an electron pair with p T > 23 GeV for the leading electron and p T > 12 GeV for the subleading electron or a muon pair with p T > 17 GeV for the leading muon and p T > 8 GeV for the subleading muon. Here leading and subleading refers to the electron or muon with the highest and second-highest p T , respectively, in the event. The different-flavour dilepton triggers require either an electron with p T > 23 GeV and a muon with p T > 8 GeV, or a muon with p T > 23 GeV and an electron with p T > 8 GeV.
Events are reconstructed using a particle-flow (PF) algorithm [62], which aims to identify and reconstruct each individual particle in an event with an optimised combination of information from the various elements of the CMS detector. Charged hadrons from pileup are subtracted on an event-by-event basis. Subsequently, the remaining neutral-hadron component from pileup is accounted for through jet energy corrections [63].
Electron candidates are reconstructed from a combination of the track momentum at the main interaction vertex, the corresponding energy deposition in the ECAL, and the energy sum of all bremsstrahlung photons associated with the track [64]. The electron candidates are required to have p T > 25 GeV for the leading candidate and p T > 20 GeV for the subleading candidate and |η| < 2.4. Electron candidates with ECAL clusters in the region between the barrel and endcap (1.44 < |η cluster | < 1.57) are excluded, because the reconstruction of an electron object in this region is not optimal. A relative isolation criterion I rel < 0.06 is applied, where I rel is defined as the p T sum of all neutral hadron, charged hadron, and photon candidates within a distance of 0.3 from the electron in η-φ space, divided by the p T of the electron candidate. In addition, electron identification requirements are applied to reject misidentified electron candidates and candidates originating from photon conversions. Muon candidates are reconstructed using the track information from the tracker and the muon system. They are required to have p T > 25 GeV for the leading candidate and p T > 20 GeV for the subleading candidate and |η| < 2.4. An isolation requirement of I rel < 0.15 is applied to muon candidates, including particles within a distance of 0.4 from the muon in η-φ space. In addition, muon identification requirements are applied to reject misidentified muon candidates and candidates originating from in-flight decay processes. For both electron and muon candidates, a correction is applied to I rel to suppress residual pileup effects.
Jets are reconstructed by clustering the PF candidates using the anti-k T clustering algorithm [65,66] with a distance parameter R = 0.4. To correct for the detector response, p T -and η-dependent jet energy corrections are applied [67]. A jet is selected if it has p T > 30 GeV and |η| < 2.4. Jets are rejected if the distance in η-φ space between the jet and the closest lepton, ∆R(jet, lepton), is less than 0.4. Jets originating from the hadronisation of b quarks (b jets) are identified with an algorithm [68] that uses secondary vertices together with track-based lifetime information to construct a b tagging discriminant. The chosen working point has a b jet tagging efficiency of ≈80-90% and a mistagging efficiency of ≈10% for jets originating from gluons, as well as u, d, or s quarks, and ≈30-40% for jets originating from c quarks.
The missing transverse momentum vector p miss T is defined as the projection on the plane perpendicular to the beams of the negative vector sum of the momenta of all PF candidates in an event. Its magnitude is referred to as p miss T . Jet energy corrections are propagated to improve the determination of p miss T .
Events are selected offline if they contain exactly two isolated electrons or muons of opposite electric charge. Furthermore, they need to contain at least two jets and at least one of these jets must be b tagged. Events with an invariant mass of the lepton pair, M( ), smaller than 20 GeV are removed in order to suppress contributions from heavy-flavour resonance decays and low-mass Drell-Yan processes. Backgrounds from Z+jets processes in the e + e − and µ + µ − channels are further suppressed by requiring |m Z − M( )| > 15 GeV, and p miss T > 40 GeV. The remaining background contribution from tW, Z+jets, W+jets, diboson and tt events from decay channels other than that of the signal are estimated from the simulation.
In this analysis, the tt production cross section is also measured as a function of the extra jet multiplicity, N jet . Extra jets (also referred to as additional jets) are jets arising primarily from hard QCD radiation and not from the top quark decays. At generator level, the extra jets are defined in dilepton tt events as jets with p T > 30 GeV, |η| < 2.4, built of particles except neutrinos using the anti-k T clustering algorithm [65,66] with a distance parameter R = 0.4, and isolated from the charged leptons (i.e. e or µ) and b quarks originating from the top quark decays by a minimal distance of 0.4 in η-φ space. The charged leptons and b quarks are taken directly after W and top quark decays, respectively. At reconstruction level the extra jets are defined in dilepton tt candidate events as jets with p T > 30 GeV, |η| < 2.4, and isolated from the leptons and b jets originating from the top quark decays by the same minimal distance in η-φ space.
The tt kinematic properties are determined from the four-momenta of the decay products using a kinematic reconstruction method [15]. The three-momenta of the neutrino (ν) and of the antineutrino (ν) are not directly measured, but they can be reconstructed by imposing the following six kinematic constraints: the conservation in the event of the total transverse momentum vector, and the masses of the W bosons, top quark, and top antiquark. The reconstructed top quark and antiquark masses are required to be 172.5 GeV. The p miss T in the event is assumed to originate solely from the two neutrinos in the top quark and antiquark decay chains. To resolve the ambiguity due to multiple algebraic solutions of the equations for the neutrino momenta, the solution with the smallest invariant mass of the tt system is taken. The reconstruction is performed 100 times, each time randomly smearing the measured energies and directions of the reconstructed leptons and jets within their resolution. This smearing procedure recovers certain events that initially yield no solution because of measurement uncertainties. The three-momenta of the two neutrinos are determined as a weighted average over all smeared solutions. For each solution, the weight is calculated based on the expected true spectrum of the invariant mass of a lepton and a b jet stemming from the decay of a top quark and taking the product of the two weights for the top quark and antiquark decay chains. All possible lepton-jet combinations in the event that satisfy the requirement on the invariant mass of the lepton and jet M b < 180 GeV are considered. Combinations are ranked based on the presence of b-tagged jets in the assignments, i.e. a combination with both leptons assigned to b-tagged jets is preferred over those with one or no b-tagged jet. Among assignments with equal number of b-tagged jets, the one with the highest sum of weights is chosen. Events with no solution after smearing are discarded. The efficiency of the kinematic reconstruction, defined as the number of events where a solution is found divided by the total number of selected tt events, is studied in data and simulation and consistent results are observed. The efficiency is about 90% for signal events. After applying the full event selection and the kinematic reconstruction of the tt system, 150 410 events are observed in the e ± µ ∓ channel, 34 890 events in the e + e − channel, and 70 346 events in the µ + µ − channel. Combining all decay channels, the estimated signal fraction in data is 80.6%. Figure 1 shows the distributions of the reconstructed top quark and tt kinematic variables and of the multiplicity of additional jets in the events. In general, the data are reasonably well described by the simulation, however some trends are visible, in particular for p T (t), where the simulation predicts a somewhat harder spectrum than that observed in data, as reported in previous differential tt cross section measurements [15, 18, 20-22, 25, 26].
The M(tt) value obtained using the full kinematic reconstruction described above is highly sensitive to the value of the top quark mass used as a kinematic constraint. Since one of the objectives of this analysis is to extract the top quark mass from the differential tt measurements, exploiting the M(tt) distribution in particular, an alternative algorithm is employed, which reconstructs the tt kinematic variables without using the top quark mass constraint. This algorithm is referred to as the "loose kinematic reconstruction". In this algorithm, the νν system is reconstructed rather than the ν and ν separately. Consequently, it can only be used to reconstruct the total tt system but not the top quark and antiquark separately. As in the full kinematic reconstruction, all possible lepton-jet combinations in the event that satisfy the requirement on the invariant mass of the lepton and jet M b < 180 GeV are considered. Combinations are ranked based on the presence of b-tagged jets in the assignments, but among combinations with equal number of b-tagged jets, the ones with the highest-p T jets are chosen. The kinematic variables of the νν system are derived as follows: its p T is set equal to p miss T , while its unknown longitudinal momentum and energy are set equal to the longitudinal momentum and energy of the lepton pair. Additional constraints are applied on the invariant mass of the neutrino pair, M(νν) ≥ 0, and on the invariant mass of the W bosons, M(W + W − ) ≥ 2M W , which have only minor effects on the performance of the reconstruction. The method yields similar tt kinematic resolutions and reconstruction efficiency as for the full kinematic reconstruction. In this analysis, the loose kinematic reconstruction is exclusively used to measure triple-differential cross sections as a function of M(tt), y(tt), and extra jet multiplicity, which are exploited to determine QCD parameters, as well as the distributions used to cross-check the results. Figure 2 shows the distributions of the reconstructed tt invariant mass and rapidity using the loose kinematic reconstruction. These distributions are similar to the ones obtained using the full kinematic reconstruction (as shown in Fig. 1).

Signal extraction and unfolding
The number of signal events in data is extracted by subtracting the expected number of background events from the observed number of events for each bin of the observables. All expected background numbers are obtained directly from the MC simulations (see Section 3) except for tt final states other than the signal. The latter are dominated by events in which one or both of the intermediate W bosons decay into τ leptons with subsequent decay into an electron or muon. These events arise from the same tt production process as the signal and thus the normalisation of this background is fixed to that of the signal. For each bin the number of  , y(t) (upper right), p T (tt) (middle left), y(tt) (middle right), M(tt) (lower left), and N jet (lower right) in selected events after the kinematic reconstruction, at detector level. The experimental data with the vertical bars corresponding to their statistical uncertainties are plotted together with distributions of simulated signal and different background processes. The hatched regions correspond to the estimated shape uncertainties in the signal and backgrounds (as detailed in Section 7). The lower panel in each plot shows the ratio of the observed data event yields to those expected in the simulation. events obtained after the subtraction of other background sources is multiplied by the ratio of the number of selected tt signal events to the total number of selected tt events (i.e. the signal and all other tt events) in simulation.
The numbers of signal events obtained after background subtraction are corrected for detector effects, using the TUNFOLD package [69]. The event yields in the e + e − , µ + µ − and e ± µ ∓ channels are added together, and the unfolding is performed. It is verified that the measurements in the separate channels yield consistent results. The response matrix plays a key role in this unfolding procedure. An element of this matrix specifies the probability for an event originating from one bin of the true distribution to be observed in a specific bin of the reconstructed observables. The response matrix includes the effects of acceptance, detector efficiencies, and resolutions. The response matrix is defined such that the true level corresponds to the full phase space (with no kinematic restrictions) for tt production at parton level. At the detector level, the same kinematic ranges are used as at the generator level, but with the total number of bins typically a few times larger. The response matrix is taken from the signal simulation. The generalised inverse of the response matrix is used to obtain the distribution of unfolded event numbers from the measured distribution by applying a χ 2 minimisation technique. An additional χ 2 term is included representing Tikhonov regularisation [70]. The regularisation reduces the effect of the statistical fluctuations present in the measured distribution on the highfrequency content of the unfolded spectrum. The regularisation strength is chosen such that the global correlation coefficient is minimal [71]. For the measurements presented here, this choice results in a small contribution from the regularisation term to the total χ 2 , on the order of a few percent. The unfolding of multidimensional distributions is performed by internally mapping the multi-dimensional arrays to one-dimensional arrays [69].

Cross section determination
The normalised cross sections for tt production are measured in the full tt kinematic phase space at parton level. The number of unfolded signal eventsM unf i in bins i of kinematic variables is used to define the normalised cross sections as a function of several (two or three) variables where the total cross section σ is evaluated by summing σ i over all bins, B is the branching ratio of tt into e + e − , µ + µ − , and e ± µ ∓ final states and L is the integrated luminosity of the data sample. For presentation purposes, the measured cross sections are divided by the bin width of the first variable. They present single-differential cross sections as a function of the first variable in different ranges of the second or second and third variables and are referred to as double-or triple-differential cross sections, respectively. The bin widths are chosen based on the resolutions of the kinematic variables, such that the purity and the stability of each bin is generally above 20%. For a given bin, the purity is defined as the fraction of events in the tt signal simulation that are generated and reconstructed in the same bin with respect to the total number of events reconstructed in that bin. To evaluate the stability, the number of events in the tt signal simulation that are generated and reconstructed in a given bin are divided by the total number of reconstructed events generated in the bin.

Systematic uncertainties
The systematic uncertainties in the measured differential cross sections are categorised into two classes: experimental uncertainties arising from imperfect modelling of the detector response, and theoretical uncertainties arising from the modelling of the signal and background processes. Each source of systematic uncertainty is assessed by changing in the simulation the corresponding efficiency, resolution, or scale by its uncertainty, using a prescription similar to the one followed in Ref. [26]. For each change made, the cross section determination is repeated, and the difference with respect to the nominal result in each bin is taken as the systematic uncertainty.

Experimental uncertainties
To account for the pileup uncertainty, the value of the total pp inelastic cross section, which is used to estimate the mean number of additional pp interactions, is varied by ±4.6%, corresponding to the uncertainty in the measurement of this cross section [60].
The efficiencies of the dilepton triggers are measured with independent triggers based on a p miss T requirement. Scale factors, defined as the ratio of the trigger efficiencies in data and simulation, are calculated in bins of lepton η and p T . They are applied to the simulation and varied within their uncertainties. The uncertainties from the modelling of lepton identification and isolation efficiencies are determined using the "tag-and-probe" method with Z+jets event samples [72,73]. The differences of these efficiencies between data and simulation in bins of η and p T are generally less than 10% for electrons, and negligible for muons. The uncertainty is estimated by varying the corresponding scale factors in the simulation within their uncertainties.
The uncertainty arising from the jet energy scale (JES) is determined by varying the twentysix sources of uncertainty in the JES in bins of p T and η and taking the quadrature sum of the effects [67]. The JES variations are also propagated to the uncertainties in p miss T . The uncertainty from the jet energy resolution (JER) is determined by the variation of the simulated JER by ± 1 standard deviation in different η regions [67]. An additional uncertainty in the calculation of p miss T is estimated by varying the energies of reconstructed particles not clustered into jets.
The uncertainty due to imperfect modelling of the b tagging efficiency is determined by varying the measured scale factor for b tagging efficiencies within its uncertainties [68].
The uncertainty in the integrated luminosity of the 2016 data sample recorded by CMS is 2.5% [74] and is applied simultaneously to the normalisation of all simulated distributions.

Theoretical uncertainties
The uncertainty arising from missing higher-order terms in the simulation of the signal process at ME level is assessed by varying the renormalisation, µ r , and factorisation, µ f , scales in the POWHEG simulation up and down by factors of two with respect to the nominal values. In the POWHEG sample, the nominal scales are defined as where p T,t denotes the p T of the top quark in the tt rest frame. In total, three variations are applied: one with the factorisation scale fixed, one with the renormalisation scale fixed, and one with both scales varied up and down coherently together. The maximum of the resulting measurement variations is taken as the final uncertainty. In the parton-shower simulation, the corresponding uncertainty is estimated by varying the scale up and down by factors of 2 for initial-state radiation and √ 2 for final-state radiation, as suggested in Ref. [48].
The uncertainty from the choice of PDF is assessed by reweighting the signal simulation according to the prescription provided for the CT14 NLO set [59]. An additional uncertainty is independently derived by varying the α S value within its uncertainty in the PDF set. The dependence of the measurement on the assumed top quark mass m MC t value is estimated by varying m MC t in the simulation by ±1 GeV around the central value of 172.5 GeV.
The uncertainty originating from the scheme used to match the ME-level calculation to the parton-shower simulation is derived by varying the h damp parameter in POWHEG in the range 0.996m MC t < h damp < 2.239m MC t , according to the tuning results from Ref. [46]. The uncertainty related to modelling of the underlying event is estimated by varying the parameters used to derive the CUETP8M2T4 tune in the default setup. The default setup in PYTHIA includes a model of colour reconnection based on MPIs with early resonance decays switched off. The analysis is repeated with three other models of colour reconnection within PYTHIA: the MPI-based scheme with early resonance decays switched on, a gluon-move scheme [75], and a QCD-inspired scheme [76]. The total uncertainty from colour reconnection modelling is estimated by taking the maximum deviation from the nominal result.
The uncertainty from the knowledge of the b quark fragmentation function is assessed by varying the Bowler-Lund function within its uncertainties [77]. In addition, the analysis is repeated using the Peterson model for b quark fragmentation [78], and the final uncertainty is determined, separately for each measurement bin, as an envelope of the variations of the normalised cross section resulting from all variations of the b quark fragmentation function. An uncertainty from the semileptonic branching ratios of b hadrons is estimated by varying them according to the world average uncertainties [43]. As tt events producing electrons or muons originating from the decay of τ leptons are considered to be background, the measured differential cross sections are sensitive to the branching ratios of τ leptons decaying into electrons or muons assumed in the simulation. Hence, an uncertainty is determined by varying the branching ratios by 1.5% [43] in the simulation.
The normalisations of all non-tt backgrounds are varied up and down by ±30% taken from measurements as explained in Ref. [73].
The total systematic uncertainty in each measurement bin is estimated by adding all the contributions described above in quadrature, separately for positive and negative cross section variations. If a systematic uncertainty results in two cross section variations of the same sign, the largest one is taken, while the opposite variation is set to zero.

Results of the measurement
The normalised differential cross sections of tt production are measured in the full phase space at parton level (after radiation and before the top quark and antiquark decays) for the following variables: 1. double-differential cross sections as a function of pair of variables: • |y(t)| and p T (t), • M(tt) and |y(t)|, • M(tt) and |y(tt)|, • M(tt) and ∆η(t, t), • M(tt) and ∆φ(t, t), • M(tt) and p T (tt) and • M(tt) and p T (t).
These cross sections are denoted in the following as [y(t), p T (t) ], etc.
The pairs of variables for the double-differential cross sections are chosen in order to obtain representative combinations that are sensitive to different aspects of the tt production dynamics, mostly following the previous measurement [20]. The variables for the triple-differential cross sections are chosen in order to enhance sensitivity to the PDFs, α S , and m pole t . In particular, the combination of y(tt) and M(tt) variables provides sensitivity for the PDFs, as demonstrated in [20], the N jet distribution for α S and M(tt) for m pole t .
The numerical values of the measured cross sections and their uncertainties are provided in Appendix A. In general, the total uncertainties for all measured cross sections are about 5-10%, but exceed 20% in some regions of phase space, such as the last N jet range of the [N 0,1,2+ jet , M(tt), y(tt) ] distribution. The total uncertainties are dominated by the systematic uncertainties receiving similar contributions from the experimental and theoretical systematic sources. The largest experimental systematic uncertainty is associated with the JES. Both the JES and signal modelling systematic uncertainties are also affected by the statistical uncertainties in the simulated samples that are used for the evaluation of these uncertainties.
In Figs. 3-11, the measured cross sections are compared to three theoretical predictions based on MC simulations: POWHEG + PYTHIA ('POW+PYT'), POWHEG + HERWIG++ ('POW+HER'), and MG5 aMC@NLO + PYTHIA ('MG5+PYT'). The 'POW+PYT' and 'POW+HER' theoretical predictions differ by the parton-shower method, hadronisation and event tune (p T -ordered parton showering, string hadronisation model and CUETP8M2T4 tune in 'POW+PYT', or angular ordered parton showering, cluster hadronisation model and EE5C tune in 'POW+HER'), while the 'POW+PYT' and 'MG5+PYT' predictions adopt different matrix elements (inclusive tt production at NLO in 'POW+PYT', or tt with up to two extra partons at NLO in 'MG5+PYT') and different methods for matching with parton shower (correcting the first parton shower emission to the NLO result in 'POW+PYT', or subtracting from the exact NLO result its parton shower approximation in 'MG5+PYT'). For each comparison, a χ 2 and the number of degrees of freedom (dof) are reported. The χ 2 value is calculated taking into account the statistical and systematic data uncertainties, while ignoring uncertainties of the predictions: where R N−1 is the column vector of the residuals calculated as the difference of the measured cross sections and the corresponding predictions obtained by discarding one of the N bins, and Cov N−1 is the (N − 1) × (N − 1) submatrix obtained from the full covariance matrix by discarding the corresponding row and column. The matrix Cov N−1 obtained in this way is invertible, while the original covariance matrix Cov is singular because for normalised cross sections one degree of freedom is lost, as can be deduced from Eq. (1). The resulting χ 2 does not depend on which bin is excluded. The covariance matrix Cov is calculated as: where Cov unf and Cov syst are the covariance matrices corresponding to the statistical uncertainties from the unfolding, and the systematic uncertainties, respectively. The systematic covariance matrix Cov syst is calculated as: where C i,k,l stands for the systematic uncertainty from variation l of source k in the ith bin, and N k is the number of variations for source k. The sums run over all sources of the systematic uncertainties and all corresponding variations. Most of the systematic uncertainty sources in this analysis consist of positive and negative variations and thus have N k = 2, whilst several model uncertainties (the model of colour reconnection and the b quark fragmentation function) consist of more than two variations, a property which is accounted for in Eq. (4). All systematic uncertainties are treated as additive, i.e. the relative uncertainties are used to scale the corresponding measured value in the construction of Cov syst . This treatment is consistent with the cross section normalisation. The cross section measurements for different multi-differential distributions are statistically and systematically correlated. No attempt is made to quantify the correlations between bins from different multi-differential distributions. Thus, quantitative comparisons between theoretical predictions and the data can only be made for each single set of multi-differential cross sections.
In Fig. 3, the p T (t) distribution is compared in different ranges of |y(t)| to predictions from 'POW+PYT', 'POW+HER', and 'MG5+PYT'. The data distribution is softer than that of the predictions over the entire y(t) range. Only 'POW+HER' describes the data well, while the other two simulations predict a harder p T (t) distribution than measured in the data over the entire y(t) range. The disagreement is strongest for 'POW+PYT'.
Figures 4 and 5 illustrate the distributions of |y(t)| and |y(tt)| in different M(tt) ranges compared to the same set of MC models. The shapes of the y(t) and y(tt) distributions are reasonably well described by all models, except for the largest M(tt) range, where all theoretical predictions are more central than the data for y(t) and less central for y(tt). The M(tt) distribution is softer in the data than in the theoretical predictions. The latter trend is the strongest for 'POW+PYT', being consistent with the disagreement for the p T (t) distribution (as shown in Fig. 3). The best agreement for both [M(tt), y(t) ] and [M(tt), y(tt) ] cross sections is provided by 'POW+HER'.
In Fig. 6, the ∆η(t, t) distribution is compared in the same M(tt) ranges to the theoretical predictions. For all generators, there is a discrepancy between the data and simulation for the medium and high M(tt) bins, where the predicted ∆η(t, t) values are too low. The disagreement is the strongest for 'MG5+PYT'.
Figures 7 and 8 illustrate the comparison of the distributions of ∆φ(t, t) and p T (tt) in the same M(tt) ranges to the theoretical predictions. Both these distributions are sensitive to gluon radiation. All MC models describe the data well within uncertainties, except for 'MG5+PYT', which predicts a p T (tt) distribution in the last M(tt) bin of the [M(tt), p T (tt) ] cross sections that is too hard.
In Fig. 9, the p T (t) distribution is compared in different M(tt) ranges to the theoretical predictions. None of the MC generators is able to describe the data, generally predicting a too hard p T (t) distribution. The discrepancy is larger at high M(tt) values where the softer p T (t) spectrum in the data must be kinematically correlated with the larger ∆η(t, t) values (as shown in Fig. 6), compared to the predictions. The disagreement is the strongest for 'POW+PYT'. While the 'POW+HER' simulation is able to reasonably describe the p T (t) distribution in the entire range of y(tt) (as shown in Fig. 3), it does not provide a good description in all ranges of M(tt), in particular predicting a too hard p T (t) distribution at high M(tt).
Figures 10 and 11 illustrate the triple-differential cross sections as a function of |y(tt)| in different M(tt) and N jet ranges, measured using two or three bins of N jet . For the [N 0,1+ jet , M(tt), y(tt) ] measurement, all MC models describe the data well. For the [N 0,1,2+ jet , M(tt), y(tt) ] measurement, only 'POW+PYT' is in satisfactory agreement with the data. In particular, 'POW+HER' predicts too high a cross section for N jet > 1, while 'MG5+PYT' provides the worst description of the M(tt) distribution for N jet = 1.
All obtained χ 2 values, ignoring theoretical uncertainties, are listed in Table 1. The corresponding p-values are visualised in Fig. 12. From these values one can conclude that none of the central predictions of the considered MC generators is able to provide predictions that correctly describe all distributions. In particular, for [M(tt), ∆η(t, t) ] and [M(tt), p T (t) ] the χ 2 values are relatively large for all MC generators. In total, the best agreement with the data is provided by 'POW+PYT' and 'POW+HER', with 'POW+PYT' better describing the measurements probing N jet and radiation, and 'POW+HER' better describing the ones involving probes of the p T distribution. To extract α S and m pole t , the measured triple-differential cross sections are compared to fixedorder NLO predictions that do not have variable parameters, except for the factorisation and renormalisation scales. These provide a more rigorous assessment of theoretical uncertainties than predictions from MC event generators that complement fixed-order computations with parton showers, because the modelling of the showers might involve different PDFs and α S values, thus complicating their interpretation. Furthermore, for PDF fits using these data (to be discussed in Section 10), fast computation techniques are required that are currently available only for fixed-order calculations.
Fixed-order theoretical calculations for fully differential cross sections for inclusive tt production are publicly available at NLO O(α 3 S ) in the fixed-flavour number scheme [79], and for tt production with one [80] and two [81,82] additional jets. These calculations are used in the present analysis. Furthermore, NLO predictions for tt production with three additional jets exist [83], but are not used in this paper because the sample of events with three additional                   jets is not large enough to allow us to measure multi-differential cross sections. The exact fully differential NNLO O(α 4 S ) calculations for inclusive tt production have recently appeared in the literature [84,85], but these predictions have not been published yet for multi-differential cross sections. The NNLO calculations for tt production with additional jets have not been performed yet. values. The NLO predictions are obtained using the MG5 aMC@NLO framework running in the fixed-order mode. The cross sections for the N jet = 0 (1) bin are obtained as the difference of the inclusive tt (tt + 1 jet) and tt + 1 jet (tt + 2 jets) cross sections. The normalisation cross section is evaluated by integrating the differential cross sections over all bins, i.e. it is given by the inclusive tt cross section. A number of the latest proton NLO PDF sets are used, namely: ABMP16 [86], CJ15 [87], CT14 [59], HERAPDF2.0 [88], JR14 [89], MMHT2014 [90], and NNPDF3.1 [91], available via the LHAPDF interface (version 6.1.5) [92]. No tt data were used in the determination of the CJ15, CT14, HERAPDF2.0 and JR14 PDF sets; only total tt production cross section measurements were used to determine the ABMP16 and MMHT2014 PDFs, and both total and differential (from LHC Run 1) tt cross sections were used in the NNPDF3.1 extraction. The number of active flavours is set to n f = 5, an m pole t = 172.5 GeV is used, and α S is set to the value used for the corresponding PDF extraction. The renormalisation and factorisation scales are chosen to be µ r = µ f = H /2, H = ∑ i m t,i . Here the sum is running over all final-state partons (t, t, and up to three light partons in the tt + 2 jet calculations) and m t denotes a transverse mass, defined as m t = √ m 2 + p 2 T . The theoretical uncertainty is estimated by varying µ r and µ f independently up and down by a factor of 2, with the additional restriction that the ratio µ r /µ f stays between 0.5 and 2 [93]. Additionally, an alternative scale choice µ r = µ f = H/2, H = ∑ i m t,i , with the sum running only over t and t [85], is considered. The scales are varied coherently in the predictions with different N jet . The final uncertainty is determined as an envelope of all scale variations on the normalised cross sections. This uncertainty is referred to hereafter as a scale uncertainty and is supposed to estimate the impact of missing higher-order terms. The PDF uncertainties are taken into account in the theoretical predictions for each PDF set. The PDF uncertainties of CJ15 [87] and CT14 [59], evaluated at 90% confidence level (CL), are rescaled to the 68% CL for consistency with other PDF sets. The uncertainties in the normalised tt cross sections originating from α S and m To compare the measured cross sections to the NLO QCD calculations, the latter are further corrected from parton to particle level. The NLO QCD calculations are provided for parton-level jets and stable top quarks, therefore the corrections (further referred to as NP) are determined using additional POWHEG + PYTHIA MC simulations for tt production with and without MPI, hadronisation and top quark decays, and defined as: Here σ particle isolated from t → , b is the cross section with MPI and hadronisation for jets built of particles excluding neutrinos and isolated from charged leptons and b quarks from the top quark decays, as defined in Section 4, and σ parton no MPI, no had., no tt decays is the cross section without MPI and hadronisation for jets built of partons excluding t and t. Both cross sections are calculated at NLO matched with parton showers. The C NP factors are used to correct the NLO predictions to particle level. The NP corrections are determined in bins of the triple-differential cross sections as a function of N jet , M(tt), and y(tt), even though they depend primarily on N jet and have only weak dependence on the tt kinematic properties. For the cross sections with up to two extra jets measured in this analysis, the estimated NP corrections are close to 1, within 5%. The dependence of the NP corrections on MC modelling was studied using MC samples with varied hadronisation model, underlying event tune, and ME and parton-shower scales, as detailed in Section 7. All resulting variations of C NP were found to be 1%, therefore no uncertainties on the determined NP corrections are assigned. To compare to the measured cross sections, the normalised multi-differential cross sections of the theoretical predictions are obtained by dividing the cross sections in specific bins by the total cross section summed over all bins.      values. For each comparison, a χ 2 is calculated, taking into account the uncertainties of the data but ignoring uncertainties of the predictions. For the comparison in Fig. 14, additional χ 2 values are determined, taking also PDF uncertainties in the predictions into account, i.e. Eq. (3) becomes Cov = Cov unf + Cov syst + Cov PDF , where Cov PDF is a covariance matrix that accounts for the PDF uncertainties. Theoretical uncertainties from scale, α S (m Z ), and m pole t variations are not included in this χ 2 calculation. Sizeable differences of the χ 2 values are observed for the predictions obtained using different PDFs. Among the PDF sets considered, the best description of the data is provided by the ABMP16 PDFs. This comparison also shows that the data prefer lower α S (m Z ) and m      Fig. 18 to the world average [94]. The contributions to the total uncertainty arising from the data and from the theoretical prediction due to PDF, scale, and m pole t or α S (m Z ) uncertainties are shown separately. For the extraction of α S (m Z ), the uncertainties from the data and from PDF, scale and m pole t are comparable in magnitude. The size of the PDF uncertainties varies significantly for different PDF sets, and the extracted α S (m Z ) values depend on the input PDFs because of a strong correlation between α S and the gluon distribution. This illustrates that precise and reliable α S (m Z ) extractions from the observed data can be obtained only in a simultaneous PDF and α S (m Z ) fit. For the m pole t extraction, the total uncertainty is dominated by the data uncertainties. The world average is computed based on extractions of m pole t from inclusive tt cross sections at NNLO+NNLL and differential distributions at NLO, and dominated by the inclusive cross section measurement and a measurement from leptonic distributions. For the combination, correlations were not taken into account.
Near the mass threshold, relevant for the m pole t extraction, the fixed-order perturbation series should be improved with all-order soft-gluon resummation that, however, is not available in the tools used to obtain theoretical predictions in this work. In Ref. [95] these effects are found to be relevant only very close to the threshold (within a few GeV) and give a correction of about +1% to the total tt cross section. Attributing a 1% correction to the first M(tt) bin and assuming it is independent of |y(tt)| and N jet , it results in a shift of +0.7 GeV for the extracted value of m pole t in this analysis, which is comparable to the total uncertainty. Furthermore, the impact of the parton shower was discussed in Ref. [7], where the predictions for tt + 1 jet produc- tion obtained at NLO and using POWHEG NLO calculations matched with the PYTHIA parton shower have been compared (as shown in Fig. 1 of Ref. [7]) and agreement between different approaches was found to be within 0.5 GeV for the extracted m pole t . In the future theoretical calculations should include resummation effects to accurately extract m pole t from differential tt cross sections. Furthermore, electroweak corrections can be significant in some regions of phase space [96,97], but for the analysis in this paper no electroweak corrections were applied.
The α S and m pole t extraction is validated by repeating the procedure: 1. Using single-differential N jet , M(tt), |y(tt)| cross sections. The plots are available in Appendix B. The largest sensitivity to α S is observed when using the |y(tt)| cross sections; the value for α S is, however, strongly dependent on the PDF set used. The N jet distribution provides a smaller α S sensitivity, but with little dependence on the PDFs. For m pole t , the largest sensitivity is observed when using the M(tt) cross sections. In fact, almost no sensitivity to m pole t is present in |y(tt)| or N jet single-differential cross sections. All determinations using the single-differential cross sections yielded α S (m Z ) and m values are consistent with the nominal ones obtained using two N jet bins and have similar precision with a slightly different uncertainty composition: smaller data uncertainties but larger scale uncertainties are present when using three N jet bins. The different uncertainties are expected since more N jet bins provide more sensitivity to α S , while the NLO theoretical prediction for the last N jet bins (two or more extra jets) have larger scale uncertainties compared to the other bins (as shown in Fig. 13). This shows that NLO QCD predictions are able to describe tt data with up to two hard extra jets, however higher-order calculations are desirable to match the experimental precision in order to achieve a most accurate α S and m pole t determination.
3. Using triple-differential [p T (tt), M(tt), y(tt)] cross sections with two p T (tt) bins. The NLO calculations for inclusive tt and tt + 1 jet production with an appropriate jet p T threshold are used to describe the distribution in the two p T (tt) bins (see Appendix B.1 for further details). The extracted α S (m Z ) and m fit at NLO (also referred to as a QCD analysis, or PDF fit), together with the combined HERA inclusive deep inelastic scattering (DIS) data [88]. The XFITTER program (version 2.0.0) [98], an open-source QCD fit framework for PDF determination, is used. The precise HERA DIS data, obtained from the combination of individual H1 and ZEUS results, are directly sensitive to the valence and sea quark distributions and probe the gluon distribution through scaling violations. Therefore, these data form the core of all PDF fits. The measured tt cross sections are included in the fit to constrain α S , m pole t and the gluon distribution at high values of x, where x is the fraction of the proton momentum carried by a parton. The typical probed x values can be estimated using the LO kinematic relation The present measurement is expected to be mostly sensitive to x values in the region 0.01 x 0.1, as estimated using the highest or lowest |y(tt)| or M(tt) bins and taking the low or high bin edge where the cross section is largest.

Details of the QCD analysis
The scale evolution of partons is calculated through DGLAP equations [99][100][101][102][103][104][105] at NLO, as implemented in the QCDNUM program [106] (version 17.01.14). The Thorne-Roberts [107][108][109] variable-flavour number scheme at NLO is used for the treatment of the heavy-quark contributions. The number of flavours is set to 5, with c and b quark mass parameters M c = 1.47 GeV and M b = 4.5 GeV [88]. For the DIS data µ r and µ f are set to Q, which denotes the fourmomentum transfer. The Q 2 range of the HERA data is restricted to Q 2 > Q 2 min = 3.5 GeV 2 [88]. The theoretical predictions for the tt cross sections are calculated as described in Section 9 and are included in the fit using the MG5 aMC@NLO (version 2.6.0) [37]  The procedure for the determination of the PDFs follows the approach of HERAPDF2.0 [88]. The parametrised PDFs are the gluon distribution xg(x), the valence quark distributions xu v (x) and xd v (x), and the u-and d-type antiquark distributions xU(x) and xD(x). At the initial QCD evolution scale µ 2 f0 = 1.9 GeV 2 , the PDFs are parametrised as: assuming the relations xU(x) = xu(x) and xD(x) = xd(x) + xs(x). Here, xu(x), xd(x), and xs(x) are the up, down, and strange antiquark distributions, respectively. The sea quark distribution is defined as and A g are determined by the QCD sum rules. The B and B parameters determine the PDFs at small x, and the C parameters describe the shape of the distributions as x → 1. The parameter C g is fixed to 25 such that the term does not contribute at large x [88,112]. Additional constraints B U = B D and A U = A D (1 − f s ) are imposed to ensure the same normalisation for the xu and xd distributions as x → 0. The strangeness fraction f s = xs/(xd + xs) is fixed to f s = 0.4 as in the HERAPDF2.0 analysis [88]. This value is consistent with the determination of the strangeness fraction when using the CMS measurements of W+c production [113].
The D and E parameters are added for some distributions in order to provide a more flexible functional form. The parameters in Eq. (7) are selected by first fitting with all D and E parameters set to zero, and then including them independently one at a time in the fit. The improvement in the χ 2 of the fit is monitored and the procedure is stopped when no further improvement is observed. This leads to a 15-parameter fit. The χ 2 definition used for the HERA DIS data follows that of Eq. (32) in Ref. [88]. It includes an additional logarithmic term that is relevant when the estimated statistical and uncorrelated systematic uncertainties in the data are rescaled during the fit [114]. For the tt data presented here, a χ 2 definition without such a logarithmic term is employed. The treatment of the experimental uncertainties in the tt double-differential cross section measurements follows the prescription given in Section 8. The correlated systematic uncertainties are treated through nuisance parameters. For each nuisance parameter a Gaussian probability density function is assumed and a corresponding penalty term is added to the χ 2 . The treatment of the experimental uncertainties for the HERA DIS data follows the prescription given in Ref. [88].
The uncertainties are estimated according to the general approach of HERAPDF2.0 [88] in which the fit, model, and parametrisation uncertainties are taken into account. Fit uncertainties are determined using the criterion of ∆χ 2 = 1. Model uncertainties arise from the variations in the values assumed for the c quark mass parameter of 1.41 ≤ M c ≤ 1.53 GeV, the strangeness fraction 0.3 ≤ f s ≤ 0.5, and the value of Q 2 min imposed on the HERA data. The latter is varied within 2.5 ≤ Q 2 min ≤ 5.0 GeV 2 , following Ref. [88]. The parametrisation uncertainty is estimated by varying the functional form in Eq. (7) of all parton distributions, with D and E parameters added or removed one at a time. Additional parametrisation uncertainties are considered by using two other functional forms in Eq. (7): with A g = 0 and E g = 0, since the χ 2 in these variants of the fit are only a few units worse than that with the nominal parametrisation. Furthermore, µ 2 f0 is changed from 1.9 to 1.6 and 2.2 GeV 2 . The parametrisation uncertainty is constructed as an envelope, built from the maximal differences between the PDFs or QCD parameters resulting from the central fit and all parametrisation variations. For the PDFs, this uncertainty is valid in the x range covered by the PDF fit to the data. The total uncertainty is obtained by adding the fit, model, and parametrisation uncertainties in quadrature. For α S and m pole t extraction, the scale uncertainties in the theoretical predictions for tt production are also considered.
Here 'fit', 'model' and 'param' denote the fit, model and parameter uncertainties discussed above. The uncertainties arising from the scale variations are estimated by repeating the fit with altered values of the scales as described in Section 9 and taking the differences with respect to the nominal result. The individual contributions to the uncertainties are listed in Table 2 where the correlation was obtained from the data uncertainties propagated to the fit. This shows that the two SM parameters can be simultaneously determined from these data to high precision with only weak correlation between them.
The global and partial χ 2 values are listed in Table 3, illustrating the consistency of the input data with the fit model. In particular, the tt data are well described in the fit. The DIS data show χ 2 /dof values slightly larger than unity, similar to what is observed and investigated in Ref. [88]. For the tt data, the full χ 2 (including uncorrelated and correlated data uncertainties) is 20 for 23 degrees of freedom. The tt cross sections are compared to the NLO predictions obtained after the fit in Fig. 19. Furthermore, in Fig. 20 the [y(t), p T (t) ] cross sections (which were not used in the fit) are compared to NLO predictions obtained using the fitted PDFs, α S and m pole t , as well as other global PDF sets. The data are in satisfactory agreement with the predictions obtained in this analysis. In particular, these predictions or predictions obtained using the ABMP16 PDF set describe the slope of p T (t) considerably better than the predictions obtained using the NNPDF3.1 PDF set, while the difference in the χ 2 values is less significant. Additionally, the predicted p T (t) slope is sensitive to the m

Parameter
Variation Fit uncertainty Total Model uncertainty Table 3: The global and partial χ 2 /dof values for all variants of the QCD analysis. The variant of the fit that uses the HERA DIS only is denoted as 'Nominal fit'. For the HERA measurements, the energy of the proton beam, E p , is listed for each data set, with the electron energy being E e = 27.5 GeV, CC and NC standing for charged and neutral current, respectively. The correlated χ 2 and the log-penalty χ 2 entries refer to the χ 2 contributions from the nuisance parameters and from the logarithmic term, respectively, as described in the text.   shallow χ 2 dependence on α S (m Z ) is present when using only the HERA DIS data, similar to the findings of the HERAPDF2.0 analysis [88]. Once the tt data are included in the fit, a distinctly sharper minimum in χ 2 is observed which coincides with the one found in the simultaneous PDF and α S (m Z ) fit given in Eq. (8).

Data sets
Both the tt and the HERA DIS data are sensitive to the α S (m Z ) value in the fit: while the constraints from the tt data seem to be dominant, the residual dependence of α S (m Z ) on the HERA DIS data may remain nonnegligible. There is no way to assess the latter quantitatively because the HERA DIS data cannot be removed from the PDF fit. However, as was investigated in the HERAPDF2.0 analysis [88], when using only HERA DIS the minima are strongly dependent on the Q 2 min threshold. As a cross-check, the extraction of α S (m Z ) was repeated for a larger threshold variation 2.5 ≤ Q 2 min ≤ 30.0 GeV 2 . In contrast to the results of Ref. [88] obtained using only HERA DIS data, when adding the tt data the extracted values of α S (m Z ) show no systematic dependence on Q 2 min and are consistent with the nominal result of Eq. (8) within the total uncertainty.
To demonstrate the added value of the tt cross sections, the QCD analysis is first performed using only the HERA DIS data. In this fit, α S (m Z ) is fixed to the value extracted from the fit using the tt data, α S (m Z ) = 0.1135, and the α S (m Z ) uncertainty of ±0.0016 is added to the fit uncertainties. Then the [N 0,1+ jet , M(tt), y(tt) ] measurement is added to the fit. The global and partial χ 2 values for the two variants of the fit are listed in Table 3.
The corresponding PDFs are compared in Fig. 22. The largest impact of the tt data is observed at x 0.1. In this region the gluon distribution lacks direct constraints in the fit using the HERA DIS data only. The impact on the valence and sea quark PDFs is expected because of the  correlations between the different distributions in the fit arising in the PDF evolution and from the momentum sum rule.
In Fig. 23 the total PDF uncertainties are shown for the two variants of the fits. A reduction of uncertainties is observed for the gluon distribution, especially at x ∼ 0.1 where the included tt data are expected to provide constraints, while the improvement at x 0.1 originates mainly from the reduced correlation between α S (m Z ) and the gluon PDF. A smaller uncertainty reduction is observed for other PDFs as well (valence and sea quark distributions), because of the correlations between the PDF distributions in the fit, as explained above. In addition to the fit uncertainty reduction, the tt data constrain the large asymmetric model uncertainty of the gluon PDF at high x. This uncertainty originates from the variation of Q 2 min in the fit, using the HERA DIS data only, because of a lack of direct constraints from these data.
In Fig. 24  values of x is shown, together with their correlations. For this plot, the asymmetric α S and m pole t uncertainties are symmetrised by taking the largest deviation, and the correlation of the fit uncertainties is assumed for the total uncertainties as well. The evolution of PDFs involves α S (m Z ), therefore PDFs always depend on the α S (m Z ) assumed during their extraction. When using only the HERA DIS data, the largest dependence on α S (m Z ) is observed for the gluon distribution. The tt data reduce this dependence, because they provide constraints on both the gluon distribution and α S , reducing their correlation. In addition, the multi-differential [N 0,1+ jet , M(tt), y(tt) ] cross sections provide constraints on m and gluon PDF can be determined simultaneously and their fitted values depend only weakly on each other. This makes future PDF fits at NNLO, once corresponding theoretical predictions for inclusive tt production with additional jets become available, very interesting.

Summary
A measurement was presented of normalised multi-differential tt production cross sections in pp collisions at √ s = 13 TeV, performed using events containing two oppositely charged leptons (electron or muon). The analysed data were recorded in 2016 with the CMS detector at the LHC, and correspond to an integrated luminosity of 35.9 fb −1 . The normalised tt cross section is measured in the full phase space as a function of different pairs of kinematic variables that describe either the top quark or the tt system. None of the central predictions of the tested Monte Carlo models is able to correctly describe all the distributions. The data exhibit softer transverse momentum p T (t) distributions than given by the theoretical predictions, as was reported in previous single-differential and double-differential tt cross section measurements. The effect of the softer p T (t) spectra in the data relative to the predictions is enhanced at larger values of the invariant mass of the tt system. The predicted p T (t) slopes are strongly sensitive to the parton distribution functions (PDFs) and the top quark pole mass m pole t value used in the calculations, and the description of the data can be improved by changing these parameters.
The measured tt cross sections as a function of the invariant mass and rapidity of the tt system, and the multiplicity of additional jets, have been incorporated into two specific fits of QCD parameters at next-to-leading order, after applying corrections for nonperturbative effects, together with the inclusive deep inelastic scattering data from HERA. When fitting only α S and m pole t to the data, using external PDFs, the two parameters are determined with high accuracy and rather weak correlation between them, however, the extracted α S values depend on the PDF set. In a simultaneous fit of α S , m pole t , and PDFs, the inclusion of the new multi-differential tt measurements has a significant impact on the extracted gluon PDF at large values of x, where x is the fraction of the proton momentum carried by a parton, and at the same time allows an accurate determination of α S and m  [10] ATLAS Collaboration, "Measurements of top quark pair relative differential cross-sections with ATLAS in pp collisions at √ s = 7 TeV", Eur. Phys. J. C 73 (2013) 2261, doi:10.1140/epjc/s10052-012-2261-1, arXiv:1207.5644.

A Measured cross sections
Tables A.1 to A.47 provide the measured cross sections, including their correlation matrices of statistical uncertainties and detailed breakdown of systematic uncertainties. The description of JES uncertainty sources can be found in Ref. [67].                                                                                           1  2  3  4  5  6  7  8  9  10  11  12  13  14  15 16     1  2  3  4  5  6  7  8  9  10  11  12  13  14  15 1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23 1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23               The NLO calculations for inclusive tt and tt + 1 jet production with an appropriate jet p T threshold are used to describe the distribution in the two p T (tt) bins. Because the final state of the NLO calculation for tt and at least one jet consists of at most two light partons, there can be up to two jets built from these partons, which balance the tt transverse momentum. Therefore e.g. for p T (tt) > 100 GeV there is at least one jet with p T > 50 GeV in the NLO calculation, and one can use the tt + 1 jet calculation requiring p T (tt) > 100 GeV without any requirement on the extra jet (if the jet p T threshold is not larger than 50 GeV). In this analysis, the jet p T > 30 GeV was found to correspond approximately to p T (tt) 50 GeV, therefore the boundary of 50 GeV was chosen to split the data into two p T (tt) bins. The minimum jet p T of 25 GeV was used in the NLO calculation for tt + 1 jet production (no selection on jet |η|), and the predicted events were required to have p T (tt) > 50 GeV. This calculation was used for the bin with p T (tt) > 50 GeV, while for the bin with p T (tt) < 50 GeV the difference of the predictions for inclusive tt production and tt + 1 jet production was used. The extracted values of α S (m Z ) and m