Measurement of the top quark polarization and $\mathrm{t\bar{t}}$ spin correlations using dilepton final states in proton-proton collisions at $\sqrt{s}=$ 13 TeV

Measurements of the top quark polarization and top quark pair ($\mathrm{t\bar{t}}$) spin correlations are presented using events containing two oppositely charged leptons (e$^+$e$^-$, e$^\pm\mu^\mp$, or $\mu^+\mu^-$) produced in proton-proton collisions at a center-of-mass energy of 13 TeV. The data were recorded by the CMS experiment at the LHC in 2016 and correspond to an integrated luminosity of 35.9 fb$^{-1}$. A set of parton-level normalized differential cross sections, sensitive to each of the independent coefficients of the spin-dependent parts of the $\mathrm{t\bar{t}}$ production density matrix, is measured for the first time at 13 TeV. The measured distributions and extracted coefficients are compared with standard model predictions from simulations at next-to-leading-order (NLO) accuracy in quantum chromodynamics (QCD), and from NLO QCD calculations including electroweak corrections. All measurements are found to be consistent with the expectations of the standard model. The normalized differential cross sections are used in fits to constrain the anomalous chromomagnetic and chromoelectric dipole moments of the top quark to -0.24 $<C_\text{tG}/\Lambda^{2}<$ 0.07 TeV$^{-2}$ and -0.33 $<C^{I}_\text{tG}/\Lambda^{2}<$ 0.20 TeV$^{-2}$, respectively, at 95% confidence level.


I. INTRODUCTION
The top quark is the heaviest known fundamental particle and has a lifetime on the order of 10 −25 s [1]. This is shorter than the quantum chromodynamic (QCD) hadronization time scale 1=Λ QCD ≈ 10 −24 s, and much shorter than the spin decorrelation time scale m t =Λ 2 QCD ≈ 10 −21 s [2] (where m t is the top quark mass). Thus, not only does the top quark decay before hadronization occurs, but also its spin information is preserved in the angular distribution of its decay products.
At the CERN LHC, top quarks are produced mostly in pairs via gluon fusion (gg → tt). The quarks are unpolarized at leading order (LO), owing to the parity-conserving nature (longitudinal polarization) and approximate time invariance (transverse polarization) of QCD interactions. In the standard model (SM), a small longitudinal polarization arises from electroweak (EW) corrections, while a small transverse polarization comes from absorptive terms at one loop (both <1% [3,4]). The spins of the top quarks and antiquarks are strongly correlated, and the configuration of spins depends on the invariant mass of the tt pair (m tt ), with like (unlike) helicity pairs dominating at low (high) m tt . This paper presents a measurement of all the independent coefficients of the top quark spin-dependent parts of the tt production density matrix, as described in Ref. [4], using events (labeled dileptonic) in which the decay of the tt pair leads to two oppositely charged leptons (e þ e − , e AE μ ∓ , or μ þ μ − ) in the final state. The analysis uses a data sample of proton-proton (pp) collision events collected by the CMS experiment at a center-of-mass (CM) energy of 13 TeV in 2016, corresponding to an integrated luminosity of 35.9 fb −1 [5]. Similar measurements have been made by the ATLAS Collaboration at ffiffi ffi s p ¼ 8 TeV [6]. Differential tt cross sections corresponding to a subset of the coefficients, and other observables sensitive to the top quark polarization and tt spin correlations, have been measured by the ATLAS and CMS Collaborations at ffiffi ffi s p ¼ 7, 8, and .
In this analysis, each coefficient is extracted from a measured normalized differential tt cross section, using the same event selection and reconstruction as described in Ref. [12]. The distributions are corrected to the parton level and extrapolated to the full phase space, using a refined unfolding procedure with no regularization bias. In addition to full statistical and systematic covariance matrices for each measured distribution, matrices are provided for the combined set of all measured bins, allowing constraints to be placed using several measured distributions simultaneously.
The absence of direct signals of beyond-the-SM (BSM) particles in the LHC data analyzed so far suggests that BSM phenomena might only be directly observed at an energy scale larger than that probed at the LHC. However, BSM physics could still indirectly manifest itself in new vertices and modified couplings. Such effects can be accommodated by adding higher-dimensional operators to the SM Lagrangian in an effective field theory (EFT) approach. The coefficients measured in this analysis are sensitive to all but one of the operators of mass dimension six relevant for hadronic tt production [4]. We set limits on contributions from these operators using simultaneous fits to the measured normalized differential cross sections, including constraints on the chromomagnetic and chromoelectric dipole moments of the top quark.

II. FORMALISM AND OBSERVABLES
The square of the matrix element for tt production and decay to two leptons (with appropriate color and spin summation implied) [13] can be written as Here, l refers to an electron or muon, R is the spin density matrix related to on-shell tt production, and ρ andρ are the decay spin density matrices for the top quark and antiquark, respectively. The narrow width of the top quark compared to its mass allows factorization of the production and decay processes.
The aim of this analysis is to study the properties of the R matrix, which is purely a function of the partonic initial state and production kinematic variables, and is therefore sensitive to BSM phenomena in tt production [13]. While the analysis is also sensitive to BSM effects in tt decays, these effects are heavily constrained [14,15], and therefore have a minimal effect on the measured distributions [4,13].
The production spin density matrix R can be decomposed in the t andt spin spaces using a Pauli matrix basis: where 1 is the 2 × 2 unit matrix, σ i are the Pauli matrices, and the first (second) matrix in each tensor product refers to the top quark (antiquark) spin space. The functionÃ determines the total tt production cross section and the top quark kinematic distributions,B AE are three-dimensional vectors of functions that characterize the degree of top quark or antiquark polarization in each direction, andC is a 3 × 3 matrix of functions that characterize the correlation between the top quark and antiquark spins.
We choose an orthonormal basis to decompose the top quark spin, where these functions have definite properties with respect to discrete symmetries [4,16]. This basis is illustrated in Fig. 1. Working in the tt CM frame, we use the helicity axisk defined by the top quark direction and the directionp of the incoming parton to define the direction perpendicular to the scattering planen ¼ ðp ×kÞ= sin Θ, where Θ is the top quark scattering angle. The direction in the scattering plane mutually perpendicular tok and the transverse axisn is given byr ¼ ðp −k cos ΘÞ= sin Θ.
We expandB AE i andC ij in terms of the orthonormal basis fk;r;ng: The coefficient functions b AE i , c ij , and c i are functions of the partonic CM energy squared s and cos Θ. They can each be classified with respect to P, CP, T, and Bose-Einstein symmetry, and their P and CP symmetry properties are summarized in Table I. The approximate CP invariance of the SM requires theC matrix to be symmetric (i.e., the CPodd coefficient functions vanish: c k ¼ c r ¼ c n ¼ 0) and the top quark and antiquark to have the same polarization coefficient functions (i.e., b þ i ¼ b − i ). The P invariance of QCD forces the P-odd coefficient functions to vanish in the absence of EW interactions, so large values are allowed only for the P-and CP-even spin correlations c kk , c rr , c nn , and c rk (and the transverse polarization coefficient functions b AE n , but these are zero at tree level in QCD by T invariance). Any deviation from these expectations would be a sign of BSM phenomena.
The Bose-Einstein symmetry of the gg initial state requires a redefinition of ther andn axes (which are odd under Bose-Einstein symmetry) to allow nonzero values of the relevant coefficient functions [4]: fk;r;ng → fk; signðcos ΘÞr; signðcos ΘÞng; i.e., we have used the sign of the cosine of the top quark scattering angle, which is odd under Bose-Einstein symmetry, to define a "forward" direction for each event.
FIG. 1. Coordinate system used for the spin measurements, illustrated in the scattering plane for Θ < π=2 (left) and Θ > π=2 (right), where the signs ofr andn are flipped at Θ ¼ π=2 as shown in Eq. (4). Thek axis is defined by the top quark direction, measured in the tt CM frame. For the basis used to define the coefficient functions in Eq. (3), the incoming particles p represent the incoming partons, while for the basis used to measure the coefficients in Eqs. (8)-(10) they represent the incoming protons.
The top quark spin cannot be measured directly, but the angular distribution of the decay products of a top quark is correlated with its spin axis [2]: where Γ is the top quark decay width, χ a is the angle between the direction of decay product a and the top quark spin axis in the top quark rest frame, and κ a is the spin analyzing power. The charged lepton has maximal spin analyzing power, κ l þ ≈ 1 [17]. For top antiquark decay, the sign is reversed: Each of the 15 coefficient functions from Eq. (3) (six b AE i and nine c ij=i 0 ) is probed by a normalized differential cross section at the parton level, using the charged lepton directions measured in the rest frames of their parent top quark and antiquark as proxies for the top quark and antiquark spins. Since the measurements are made in pp collisions, the basis is adjusted from that of Eq. (4) by definingp ¼ ð0; 0; 1Þ, the direction of the proton beam in the positive z direction in the laboratory frame, in the derivation ofr andn [4]. This basis is illustrated in Fig. 1. The fourfold angular distribution for the two leptons follows from Eqs. (1) and (2), and is given by where σ is the tt production cross section, Ω 1;2 are the solid angles of the leptons in their parent top quark and antiquark rest frames, andl 1;2 are the corresponding unit vectors. The negative sign in front of the matrix C is chosen to define same-helicity top quarks as having positive spin correlation. The elements of the vectors B 1;2 and the matrix C are the following coefficients [in whose definitions the factors of κ l þ and κ l − from Eq. (5) are absorbed]: (1) B i 1 and B i 2 , the top quark and antiquark polarization coefficients with respect to each reference axis i (sensitive to b þ i and b − i ).
(2) C ii , the "diagonal" spin correlation coefficient for each reference axis i (sensitive to c ii ). (3) C ij , the "cross" spin correlation coefficients for each pair of axes i ≠ j, whose sums and differences C ij AE C ji are sensitive to c ij and c i 0 . These measurable coefficients are closely related to the production spin density matrix coefficient functions from Eq. (3), but are not identical, owing to the different basis used for the spin measurement. We do not measure the coefficients differentially or attempt to separate the contributions from different initial states. The association TABLE I. Observables and their corresponding measured coefficients, production spin density matrix coefficient functions, and P and CP symmetry properties. For the laboratory-frame asymmetries shown in the last two rows, there is no direct correspondence with the coefficient functions.

Observable
Measured coefficient Coefficient function Symmetries C kk c kk P-even, CP-even cos θ r 1 cos θ r 2 C rr c rr P-even, CP-even cos θ n 1 cos θ n 2 C nn c nn P-even, CP-even cos θ r 1 cos θ k 2 þ cos θ k 1 cos θ r 2 C rk þ C kr c rk P-even, CP-even cos θ r 1 cos θ k 2 − cos θ k 1 cos θ r 2 C rk − C kr c n P-even, CP-odd cos θ n 1 cos θ r 2 þ cos θ r 1 cos θ n 2 C nr þ C rn c nr P-odd, CP-even cos θ n 1 cos θ r 2 − cos θ r 1 cos θ n 2 C nr − C rn c k P-odd, CP-odd cos θ n 1 cos θ k 2 þ cos θ k 1 cos θ n between the measured coefficients and the coefficient functions is given in Table I. For each of the 15 coefficients that make up B 1;2 and C, a change of variables can be made to obtain a singledifferential cross section that depends only on that coefficient. After integrating out the azimuthal angles, Eq. (6) reduces to where θ i 1 (θ j 2 ) is the angle of the positively (negatively) charged lepton, measured with respect to axis i (j) in the rest frame of its parent top quark (antiquark). By changing variables (if necessary) and integrating out one of the angles, we can derive single-differential cross sections with respect to cos θ i 1 , cos θ i 2 , and cos θ i 1 cos θ j 2 : Similarly, starting from Eq. (6), we obtain (for i ≠ j) Thus, in order to determine the 15 coefficients, the tt production cross section is measured as a function of each of the following 15 observables at the parton level: (1) The three cos θ i 1 terms and three cos θ i 2 terms to measure B i 1 and B i 2 , the top quark and antiquark polarization coefficients with respect to each reference axis i.
(2) The three cos θ i 1 cos θ i 2 terms to measure C ii , the diagonal spin correlation coefficient for each axis i.
(3) The six sum and difference terms cos θ i 1 cos θ j 2 AE cos θ j 1 cos θ i 2 to measure C ij AE C ji , the sums and differences of the cross spin correlation coefficients for each pair of axes i ≠ j. We do not measure the separate cross spin correlation cos θ i 1 cos θ j 2 distributions, because it is the sums and differences that are sensitive to the c ij and c i coefficients of Eq. (3) (see Table I).
In addition, we measure four further cos θ i 1;2 distributions based on modified axesk Ã andr Ã , equal to AEk or AEr, depending on the sign of jy t j − jy t j, the difference of the moduli of the top quark and antiquark rapidities in the laboratory frame. The use of the modified axes probes the coefficient functions in different areas of phase space, providing sensitivity to different combinations of fourquark operators [4].
The spin correlation coefficient D is related to the diagonal C coefficients as D We make a direct measurement of the D coefficient using the distribution of the dot product of the two lepton directions measured in their parent top quark and antiquark rest frames [4],l 1 ·l 2 ¼ cos φ: We also measure two related laboratory-frame distributions, using the following observables: (1) cos φ lab ¼l lab 1 ·l lab 2 , defined by analogy to cos φ, but using the lepton directions measured in the laboratory frame, which have excellent experimental resolution.
(2) jΔϕ ll j, the absolute value of the difference in azimuthal angle ϕ between the two leptons in the laboratory frame. The association between the 22 measured observables and the coefficients and coefficient functions is given in Table I. Except for the laboratory-frame distributions, the distribution shapes are completely determined by the coefficients, following the functional forms of Eqs. (8)-(10). The laboratory-frame observables (given in the last two rows of Table I) do not directly relate to any of the coefficients, and we instead quantify the shapes of their distributions by calculating the asymmetry (A) in the number of events (N) about the center of the distribution:

III. 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 end cap sections. Forward calorimeters extend the pseudorapidity (η) coverage provided by the barrel and end cap detectors. Muons are measured in gas-ionization detectors embedded in the steel flux-return yoke outside the solenoid. Events of interest are selected using a two-tiered trigger system [18]. 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 optimized 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. [19].

IV. EVENT SIMULATION
Simulated tt events with a top quark mass of m t ¼ 172.5 GeV are produced at next-to-leading order (NLO) in QCD at the matrix element (ME) level by the POWHEG v. 2 [20][21][22][23] generator (POWHEGv2). The h damp parameter of POWHEGv2, which regulates the damping of real emissions in the NLO calculation when matching to the parton shower (PS), is set to 272.72 GeV [24]. The PS and hadronization are performed by PYTHIA 8.219 [25] (referred to as PYTHIA8 in the following) with the CUETP8M2T4 tune [24,26,27]. In order to assess the level of variation when using an alternative ME and matching procedure, an alternative tt sample is generated using the MADGRAPH5_aMC@NLO 2.3.3 [28] generator including up to two extra partons at the ME level with NLO precision. The decays of the top quarks are modeled using MADSPIN [29], and events are matched to PYTHIA8 for PS and hadronization using the FxFx jet merging prescription [30]. This sample is referred to as MG5_aMC@NLO þ PYTHIA8 [FxFx].
Signal tt events are defined as those with two charged leptons (e þ e − , e AE μ ∓ , or μ þ μ − ), originating from W boson decays and not from τ lepton decays. All other tt events are regarded as a background. The largest background contributions originate from tt events with leptonically decaying τ leptons, single top quarks produced in association with a W boson (tW), and, in events with same-flavor leptons, Z=γ Ã bosons produced with additional jets (Z þ jets). Additional significant backgrounds include W boson production with additional jets (W þ jets), diboson (WW, WZ, and ZZ) events, and the production of a tt pair in association with a W or a Z boson (tt þ W=Z). Other sources of background are negligible in comparison to the uncertainties in the main backgrounds, and are not included in this analysis.
The W þ jets process is simulated at LO precision using MADGRAPH5_aMC@NLO with up to four additional partons at the ME level and matched to PYTHIA8 using the MLM jet merging prescription [31]. The Z þ jets process is simulated at NLO precision using MADGRAPH5_aMC@NLO with up to two additional partons at the ME level and matched to PYTHIA8 using the FxFx prescription. The tt þ W=Z processes are simulated with MADGRAPH5_aMC@NLO with NLO precision at the ME level and matched to PYTHIA8. In the case of tt þ W, one extra parton is simulated at the ME level and the calculation is matched to PYTHIA8 using the FxFx prescription. Single top quark production is simulated with POWHEG v. 1 [32,33] with the h damp parameter set to 172.5 GeV and using the CUETP8M2T4 tune in PYTHIA8. Diboson events are simulated at LO with PYTHIA8. The NNPDF3.0_lo_as_0130 and NNPDF3.0_ nlo_as_0118 [34,35] parton distribution function (PDF) sets are used for the LO and NLO simulations, respectively.
The cross sections used to normalize the simulated predictions are calculated at the highest orders of perturbative QCD currently available: next-to-next-to-leading order (NNLO) for W þ jets and Z þ jets [36]; approximate NNLO for single top quark production in the tW channel [37]; and NLO for diboson [38] and tt þ W=Z [39]. The tt simulation is normalized to a cross section of 831.8 þ19.8 −29.2 ðscaleÞ AE 35.1ðPDF þ α S Þ pb (where α S is the strong coupling constant), calculated with the TOP++2.0 program [40] at NNLO, including resummation of nextto-next-to-leading-logarithmic soft-gluon terms and assuming m t ¼ 172.5 GeV.
Additional pp interactions within the same or nearby bunch crossings (pileup) are simulated for all samples, using a pileup multiplicity distribution that reflects the distribution of reconstructed vertices in data. The interactions of particles with the CMS detector are simulated using GEANT4 (v. 9.4) [41].

V. EVENT SELECTION
The event selection [12] targets the dileptonic decay tt → bl þ νbl −ν . To maximize the trigger efficiency, both single-lepton and dilepton trigger paths are used. For the single-electron (single-muon) trigger, a transverse momentum threshold of p T ¼ 27ð24Þ GeV is applied. The sameflavor dilepton triggers require either an electron pair with p T > 23ð12Þ GeV for the leading (trailing) electron, or a muon pair with p T > 17ð8Þ GeV for the leading (trailing) muon, where leading (trailing) refers to the electron or muon with the highest (second-highest) p T in the event. The different-flavor dilepton triggers require either an electron with p T > 12 GeV and a muon with p T > 23 GeV, or an electron with p T > 23 GeV and a muon with p T > 8 GeV.
The events selected by the HLT are reconstructed offline using a particle-flow algorithm [42], which aims at reconstructing each individual particle in an event using an optimized combination of information from the various elements of the CMS detector. Electron candidates are reconstructed from a combination of the track momentum at the main interaction vertex and the corresponding clusters in the ECAL with a Gaussian sum filter algorithm [43]. Electron candidates with ECAL clusters in the region between the barrel and end cap (1.44 < jη cluster j < 1.57) have a reduced reconstruction efficiency and are excluded. A relative isolation criterion I rel < 0.0588 (0.0571) is applied for electron candidates in the barrel (end cap). The 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 candidate in η-ϕ space, divided by the p T of the electron candidate, with a correction to suppress the residual effect of pileup. Additional electron identification requirements are applied to reject misidentified electron candidates and candidates originating from photon conversions [42,43]. Muon candidates are reconstructed using the track information from the tracker and the muon system [44]. A relative isolation requirement of I rel < 0.15 within a distance of 0.4 in η-ϕ space from the muon candidate is applied. In addition, muon identification requirements are used to reject misidentified muon candidates and candidates originating from decay-in-flight processes [44]. Both electron and muon candidates are required to have p T > 25ð20Þ GeV for the leading (trailing) candidate and jηj < 2.4.
Jets are reconstructed by clustering the particle-flow candidates using the anti-k T clustering algorithm with a distance parameter of 0.4 [45,46]. The jet momentum is determined as the vector sum of all particle momenta in the jet, and is found from simulation to be within 5-10% of the true momentum over the whole p T spectrum and detector acceptance. Pileup can contribute additional tracks and calorimetric energy depositions to the jet momentum. To mitigate this effect, tracks identified to originate from pileup vertices are discarded, and an offset correction is applied to correct for remaining contributions from neutral particles from pileup [42]. Jet energy corrections are derived from simulation to bring the measured response of jets to that of particle-level jets on average. In situ measurements of the momentum imbalance in dijet, photon þ jet, Z þ jet, and multijet events are used to account for any residual differences in the jet energy scale (JES) between data and simulation. Additional selection criteria are applied to remove badly reconstructed jets [42,47]. Jets are selected if they have p T > 30 GeV and jηj < 2.4. Jets are rejected if the distance in η-ϕ space between the jet and the closest lepton is less than 0.4. Jets originating from the hadronization of b quarks (b jets) are identified (b tagged) by combining information related to secondary decay vertices reconstructed within the jets and track-based lifetime information in an algorithm (CSVv2) that provides a b jet identification efficiency of 79-87% and a probability to misidentify light-and charm-flavor jets as b jets of approximately 10 and 40%, respectively [48].
The missing transverse momentum vector ⃗p miss T is defined as the projection on the plane perpendicular to the beam axis of the negative vector sum of the momenta of all reconstructed particles in an event. Its magnitude is referred to as p miss T .
The selected events are required to have exactly two isolated electrons or muons of opposite electric charge and at least two jets. At least one of the jets is required to be b tagged. Events with a lepton-pair invariant mass m ll < 20 GeV are removed in order to suppress contributions from heavy-flavor resonance decays and low-mass Drell-Yan processes. In the e þ e − and μ þ μ − channels, backgrounds from Z þ jets processes are further suppressed by requiring p miss T > 40 GeV and vetoing events with 76 < m ll < 106 GeV. The remaining background yield from Z þ jets events, which is large in the e þ e − and μ þ μ − channels, is determined by applying a factor derived from simulation to the number of Z þ jets events observed in data in a control region where m ll is close to the Z boson mass [49,50]. A correction to account for non-Z þ jets backgrounds in the control region is derived from the e AE μ ∓ channel. The simulated Z þ jets yield is corrected by up to 5% in each channel to match the determination from data.
The four-momenta of the top quark and antiquark in each event are estimated using a kinematic reconstruction algorithm [12,49]. The algorithm considers all possible assignments of reconstructed jets and leptons to the b quarks and leptons from top quark decay, and solves for the unknown neutrino momenta using the following assumptions and constraints: p miss T is assumed to originate solely from the two neutrinos; the invariant mass of each reconstructed W boson (m W ) must equal 80.4 GeV [1]; and the invariant mass of each reconstructed top quark must equal 172.5 GeV. Effects of detector resolution are accounted for by randomly smearing the measured energies and directions of the reconstructed jets and leptons according to their simulated resolutions. The assumed value of m W is varied according to a simulated Breit-Wigner distribution, with a width of 2.1 GeV [1]. For a given application of the smearing, the solution of the equations for the neutrino momenta yielding the smallest reconstructed m tt is chosen. For each solution, a weight is calculated based on the spectrum of the true invariant mass of the lepton and b jet system from top quark decay at the particle level [12]. The weights are summed over 100 applications of the smearing, and the top quark and antiquark four-momenta are calculated as a weighted average. Considering only the combinations with the most b-tagged jets, the assignment of jets and leptons that yields the maximum sum of weights is chosen. The efficiency of the kinematic reconstruction, defined as the fraction of the selected tt events where a solution is found, is about 90% in both data and simulation. Events with no real solution for the neutrino momenta are excluded from further analysis.
After applying the full event selection, 34 890 events in the e þ e − channel, 70 346 events in the μ þ μ − channel, and 150 410 events in the e AE μ ∓ channel are observed. The difference in the e þ e − and μ þ μ − channel yields is attributable to the lower efficiencies of the electron identification and isolation requirements. The differential cross section measurements are made using the combination of events from the three channels, where the fraction of signal events in the data sample, estimated from simulation, is 79%.
In Figs. 2 and 3, distributions of all the reconstructed angular observables (defined in Sec. II) are shown. There is reasonable agreement between the data and the sum of the expected signal and background contributions given the systematic uncertainties. In addition to the systematic uncertainties discussed in Sec. VII, two uncertainties that affect only the normalization of the measured differential cross section are considered: the 2.5% uncertainty in the integrated luminosity of the data sample [5] is applied to the normalization of all simulated predictions, and a 1.5% normalization uncertainty is applied to the tt prediction to account for the uncertainty in the dileptonic branching fraction (BF) [1]. The shapes of the reconstructed distributions differ substantially from the expected parton-level functional forms of Eqs. (8)-(10) owing to the effects of bin migration, detector acceptance and efficiency, and background events, which are described in Sec. VI.

VI. UNFOLDING THE DIFFERENTIAL CROSS SECTIONS
The effects of detector acceptance and efficiency sculpt the reconstructed distributions, and the smearing introduced by the detector response, kinematic reconstruction algorithm, PS, and hadronization leads to the migration of events across bins. In order to measure the differential cross sections at the parton level in the full phase space, these effects are accounted for by using the TUnfold regularized unfolding method [51]. The response matrix used in the unfolding is calculated for each measured distribution using the default tt simulation, where the momenta of the partonlevel top quarks are defined after QCD radiation has been simulated but before the top quark decays.
To keep the bin-to-bin migrations small (to avoid strong bin-to-bin correlations in the unfolded distributions), the widths of the measurement bins are chosen according to the reconstruction resolution of the observable. This is quantified both directly by comparing the generator level and detector level observables in simulation, and by measuring the purity and stability. Purity is defined as the fraction of events in a given bin at the detector level that originate from the same bin at the generator level, and stability is defined as the fraction of events in a given bin at the generator level that are reconstructed in the same bin at the detector level. For all observables measured in the top quark rest frame, the use of six bins of uniform width is found to be well matched to the reconstruction resolution. The purities and stabilities are typically 40%. For the observables measured in the laboratory frame (cos φ lab and jΔϕ ll j), six uniformwidth bins are also used. These observables have excellent experimental resolution, and the purities and stabilities are >99%.
The presence of background events is accounted for prior to performing the unfolding. After subtracting all other background components, the background from dileptonic tt events with leptonically decaying τ leptons is subtracted as a fraction of the total remaining events. The fraction is evaluated per bin as the ratio of the background to the total dileptonic tt events in simulation. Thus, the shapes of the distributions for dileptonic tt events are taken from data, and any dependence on the total cross section used in the normalization of the simulated tt sample is avoided.
In TUnfold, a procedure based on matrix inversion is used to obtain an unfolded distribution from the measured distribution by applying a χ 2 minimization technique. The potential large statistical fluctuations and strong anticorrelations between adjacent bins arising from the matrix inversion are suppressed by introducing a term in the χ 2 expression that smooths (regularizes) the shape of the unfolded distribution [51]. The regularization term penalizes the curvature of a vector constructed from the product of the difference between the unfolded and simulated bin values and a factor calculated using the expected functional form [Eqs. (8)-(10)] such that a deviation in the coefficient corresponds to a linear change in the vector. Since linear changes are unconstrained by regularization of the curvature, and the functional forms at the parton level (which are unaffected by BSM phenomena in tt production) depend only on the coefficient, this ensures that the regularization cannot introduce a bias in the unfolded distribution. For the laboratory-frame distributions there are no such simple functional forms, and no factor is applied to the difference vector. However, this choice is of little consequence because the regularization is very weak owing to the low level of bin migration.
The use of wide bins for the response matrix loses information about its dependence inside each bin, meaning the unfolding can be biased if the physical process density differs from the simulation. Since the curvature regularization is unbiased, we make use of narrower bins in the TUnfold χ 2 minimization; a factor 4 narrower is found to be sufficient to reduce the bias from binning to a negligible level. We have thus replaced the biased implicit regularization from binning with an unbiased regularization of the curvature within each of the original bins.
The regularization level is determined for each distribution by minimizing the average global correlation coefficient (ρ avg ) [51], where ρ avg is determined after rebinning to the original six bins.
For each measured bin, we perform tests using pseudodata to confirm a linear response of the method to variations in the coefficient, and confirm that the distribution of the difference between the nominal bin value and that measured in pseudodata, normalized to the measured uncertainty, is consistent with having zero mean and unit width.
The data in the three channels are combined before unfolding in order to model correlations between channels  in situ and reduce statistical uncertainties in poorly populated regions of the response matrix. After unfolding, each distribution is normalized to unit area to measure the normalized differential cross section.

VII. SYSTEMATIC UNCERTAINTIES
The systematic uncertainties arising from the detector performance and the modeling of the signal and background processes are evaluated from the difference between the nominal measurement and that obtained by repeating the unfolding procedure using simulated events with the appropriate systematic variation. Each source of systematic uncertainty is represented by a covariance matrix for the bins of the measured normalized differential cross sections. The total systematic uncertainty is derived from the sum of these covariance matrices. In this section, each of the applied variations is detailed and categorized into experimental and theoretical sources of uncertainty.

A. Experimental sources of uncertainty
Many of the experimental sources of uncertainty relate to the scale factors (SFs), defined as the ratio of the efficiencies in data and simulation, that are applied to the simulation in order to accurately model the data.
The efficiencies of the triggers in data are measured as the fraction of events passing alternative triggers based on a p miss T requirement that also satisfy the criteria of the trigger of interest [12,52]. As the efficiency of the p miss T requirement is only weakly correlated with the dilepton trigger efficiencies, the bias introduced by the p miss T requirement is negligible. The efficiencies are close to unity in both data and simulation, as are the corresponding SFs. To estimate the uncertainty from the modeling of the trigger efficiency, the SFs are varied within their uncertainties, both globally for all bins and depending on the η of the leptons. The total trigger uncertainty is derived by taking the maximum deviation produced by the two variations in each unfolded bin.
The SFs for the lepton identification and isolation efficiencies are determined with a tag-and-probe method using Z þ jets event samples [50,53]. Measured in bins of η and p T , the SFs are generally within 10% of unity for electrons, and consistent with unity for muons. The lepton identification and isolation uncertainty is estimated by varying the SFs within their uncertainties. The efficiency of the kinematic reconstruction of the top quarks is found to be consistent between data and simulation within around 0.2%. An associated uncertainty is derived by varying the corresponding SFs by AE0.2%.
The uncertainty from the modeling of the number of pileup events is obtained by changing the inelastic pp cross section assumed in simulation by AE4.6%, consistent with the cross section uncertainty presented in Ref. [54].
The uncertainty arising from the imperfect modeling of the b tagging efficiency is determined by varying the measured SFs within their uncertainties, both globally and depending on the p T and η of the b jets. The total uncertainty is derived by taking the maximum observed deviation in each unfolded bin. The b tagging uncertainties for heavy-flavor (b and c) and light-flavor (u, d, s, and gluon) jets are calculated separately, and combined in quadrature to give the total b tagging uncertainty. To avoid double counting of the uncertainty related to the b tagging efficiency, when necessary, the SFs for b tagging efficiency are recalculated in the evaluation of the remaining sources of experimental and theoretical uncertainty described in this section, using the procedure given in Ref. [12].
The uncertainty arising from the JES is determined by varying the individual sources of uncertainty in the JES in bins of the jet p T and η, and taking the quadrature sum of the differences [55]. The JES variations are propagated to the uncertainties in p miss T . An additional uncertainty in the calculation of p miss T is estimated by varying the energies of reconstructed particles not clustered into jets (unclustered energy). The uncertainty from the jet energy resolution (JER) is determined by the variation of the JER in simulation within its uncertainty in different η regions [55].

B. Theoretical sources of uncertainty
The uncertainty arising from the missing higher-order terms in the simulation of the signal process at the ME level is assessed by varying the renormalization and factorization scales (μ R and μ F ) in the POWHEGv2 simulation up and down by a factor of 2 with respect to their nominal values, both individually and simultaneously (six variations in total).

The nominal choice for the scales is
where p T;t denotes the p T of the top quark in the tt rest frame. In the PS simulation, the corresponding uncertainty is estimated by four additional variations: changing the scale of initial-and final-state radiation individually up and down by factors of 2 and ffiffi ffi 2 p , respectively, as suggested in Ref. [27]. The total scale uncertainty is taken as the maximum deviation from the nominal prediction from all ten variations.
The uncertainty originating from the scheme used to match the ME-level calculation to the PS simulation is derived by varying the h damp parameter in POWHEGv2 by factors of 1.42 and 0.63, according to the results of tuning this parameter from Ref. [24].
The default setup in PYTHIA8 includes a multiple parton interaction (MPI) based model of color reconnection (CR) with early resonance decays switched off. To estimate the uncertainty from this choice of model, the analysis is repeated with three other CR models within PYTHIA8: the MPI-based scheme with early resonance decays switched on, a gluon-move scheme [56], and a QCD-inspired scheme [57]. The total uncertainty from CR modeling is estimated by taking the maximum deviation from the nominal result. The uncertainty related to modeling of the underlying event is estimated by varying the parameters used to derive the CUETP8M2T4 tune in the default setup.
The uncertainty from the b quark fragmentation function is assessed from the largest deviation when varying the Bowler-Lund function within its uncertainties [58] and repeating the analysis with the Peterson model for b quark fragmentation [59]. An uncertainty from the semileptonic BF of b hadrons is estimated by correcting the tt simulation to match the BF in Ref. [1].
The uncertainty from the PDFs is assessed from the standard deviation of the result when using the replicas of the NNPDF3.0 PDF set in the signal simulation [34,60]. An additional uncertainty is derived by varying the α S value within its uncertainty in the PDF set [60]. The dependence of the measurement on the assumed m t value is estimated by varying the chosen m t in the default setup by AE1 GeV with respect to the default value of 172.5 GeV.
Previous CMS studies have shown that the p T distribution of the top quark measured from data is softer than that in the NLO simulation of tt production [12,49,[61][62][63]. This is understood to arise at least partly from the missing higher-order QCD terms [64][65][66][67]. The change in the measurement when reweighting the simulated tt event sample to match the top quark p T spectrum in data is taken as a two-sided systematic uncertainty associated with the signal modeling.
Since tt events producing electrons or muons originating from the decay of τ leptons are considered a background, the measured differential cross sections are sensitive to the relative BFs of W bosons decaying to τ leptons and electrons or muons, and the τ semileptonic BFs assumed in the simulation. An uncertainty of 2.5% is assigned to the relative normalization of this background process. The shape and absolute normalization of this process is taken from data, as described in Sec. VI. The normalizations of all other backgrounds are varied by AE30% [12,50].

A. Normalized differential cross sections
Normalized differential cross sections at the parton level are measured for the 22 observables introduced in Sec. II. For the top quark polarization observables measured using the nominal (k,r andn) and modified (k Ã andr Ã ) reference axes, the results are shown in Figs. 4 and 5, respectively. The results for the observables that probe the diagonal tt spin correlation coefficients and for the laboratory-frame spin correlation observables are shown in Fig. 6. For the cross spin correlation observables, the results are shown in Fig. 7. The measured distributions are compared with predictions from the POWHEGv2 and MADGRAPH5_aMC@NLO simulations and with calculations for tt production at NLO in QCD with EW corrections [3,4], as well as similar calculations in the absence of top quark polarization or spin correlations. For the observables measured in the top quark rest frame, the latter are equivalent to the predictions of Eqs. (8)-(10), with the coefficients set to zero. For the laboratory-frame observables, dedicated calculations were made using the computational setup described in Refs. [3,68]. In addition, the only NNLO QCD prediction [69] is shown for the jΔϕ ll j distribution in Fig. 6.
The effect of spin correlations is most clearly visible in the cos φ distribution in Fig. 6, where the data strongly favor the predictions with spin correlations compared to the uncorrelated prediction. The presence of spin correlations can also be seen in the other distributions sensitive to Pand CP-even spin correlations: the three cos θ i 1 cos θ i 2 distributions and the two laboratory-frame distributions (cos φ lab and jΔϕ ll j) in Fig. 6, and the cos θ r 1 cos θ k 2 þ cos θ k 1 cos θ r 2 distribution in Fig. 7. However, the measurements are not sensitive to the small level of top quark polarization predicted in the SM, and do not significantly disfavor the unpolarized predictions in Figs. 4 and 5.
The statistical and systematic correlation matrices for the normalized differential cross sections are determined simultaneously for all 132 measured bins to allow the fitting of multiple distributions, and are shown in Fig. 8. The statistical correlations are estimated using a bootstrap resampling of the data [70], and the systematic correlations are estimated by simultaneously evaluating the systematic variations described in Sec. VII for all measured bins. The statistical correlations among bins from the same distribution exhibit a typical pattern of correlation and anticorrelation arising from the unfolding. The statistical correlations between bins from different distributions are typically small, but the relationships between some of the distributions result in stronger correlations (for example, cos θ iÃ 1;2 and cos θ i 1;2 are the same up to a sign). The systematic correlations are in general much stronger, and the pattern of positive and negative correlations reflects the relative changes in shape of the different distributions in response to the systematic variations.
The agreement between the measured distributions and the four theoretical predictions shown in Figs. 4-7 is quantified by evaluating the χ 2 , taking the uncertainties from the sum of the measured statistical and systematic covariance matrices (and not including any uncertainties in the prediction). The results are shown in Table II. For the observables measured in the top quark rest frame, there is generally good agreement between the measured distributions and all the predictions in the presence of spin correlations. For the two observables measured in the laboratory frame, there is greater variation between the predictions. The POWHEGv2 prediction best describes the data for cos φ lab , while for jΔϕ ll j the MADGRAPH5_aMC@NLO prediction provides the best agreement. The NNLO QCD prediction shown in Fig. 6 also describes the observed jΔϕ ll j distribution well, with a

B. Coefficients
From each measured normalized differential cross section, we extract the corresponding coefficient, using the functional forms of Eqs. (8)-(10) and combining the information from the measured bins in a way that minimizes the uncertainty in the coefficient. For the laboratory-frame observables, the shapes of the distributions are instead quantified by the asymmetries defined in Eq. (11). The results for all quantities are shown with their total uncertainties in Table III, where they are compared with predictions from POWHEGv2 and MADGRAPH5_aMC@NLO simulations and the NLO calculations [3,4]. The uncertainties in the NLO calculations come from varying μ R and μ F simultaneously up and down by a factor of 2 from their nominal value of m t . The NNLO QCD prediction for A jΔϕ ll j is 0.115 þ0.005 −0.001 [69], where the uncertainties are taken from the largest deviations when varying μ R and μ F individually and simultaneously up and down by a factor of 2 from the nominal choice of ðm t T þ m¯t T Þ=4. The results are also shown in Figs. 9-11. There is good agreement between the measured coefficients and all the SM predictions, while substantial variation is seen in the predicted laboratory-frame asymmetries, which have sizable scale uncertainties. In the fixed-order calculations of Refs. [3,4], the numerator and denominator of the normalized differential cross section are computed at NLO QCD including EW corrections, and then the ratio is consistently expanded to NLO. On the   FIG. 6. Unfolded data (points) and predicted (horizontal lines) normalized differential cross sections for the diagonal spin correlation observables (first two rows) and the laboratory-frame observables (bottom row). The vertical lines on the points represent the total uncertainties, with the statistical components indicated by horizontal bars. The ratios of various predictions to the data are shown in the lower panels.   (right) correlation matrices for all measured bins of the normalized differential cross sections. Each group of six bins along each axis corresponds to a measured distribution, and for conciseness is labeled by the name of the associated coefficient (as defined in Table I).
TABLE II. The χ 2 between the data and the predictions for all measured normalized differential cross sections (Figs. 4-7). The last column refers to the prediction in the case of no spin correlation or polarization. The χ 2 values are evaluated using the sum of the measured statistical and systematic covariance matrices. The number of degrees of freedom (d.o.f.) is 5 for all observables. In the last row, the χ 2 values are given for the set of all measured bins.  other hand, in the computations from simulation the ratio is not expanded. This leads to differences that are nominally of order α S 2 , in addition to the EW corrections (which are not included in the simulation). However, for A jΔϕ ll j the EW corrections are found to be only at the level of 2% using the computational setup of Ref. [3]. In the NNLO QCD calculation, the ratio is not expanded [69].
The breakdown of the systematic and statistical uncertainties in the polarization and spin correlation measurements is given in Tables IV and V, respectively. The systematic and statistical uncertainties are of comparable size for most of the measured coefficients. The exception is B k 1;2 , because the JES and b quark fragmentation uncertainties have a large effect on the reconstructed top quark momentum in the tt CM frame. The laboratory-frame asymmetries have statistical uncertainties smaller than their systematic uncertainties. The excellent reconstruction resolution results in little dilution of the statistical precision of the measured asymmetries. There is also a large background uncertainty in A lab cos φ , owing to the large Z þ jets contribution near cos φ lab ¼ 1, and a large top quark p T modeling uncertainty in A jΔϕ ll j . The statistical and systematic correlation matrices for the measured coefficients are shown in Fig. 12. The coefficients are largely statistically uncorrelated, as expected for the measurement of independent quantities. The expected statistical correlations between the related D and diagonal C coefficients are clear, as are the correlations between D and the two related laboratory-frame asymmetries. The systematic correlations are in general much stronger. In particular, strong correlations are evident for the polarization measurements with positively and negatively charged leptons (except for the B n i , where the largest sources of systematic uncertainty have a substantial statistical uncertainty from the simulation). The coefficients with significant statistical correlations naturally have significant systematic correlations as well.
The sums and differences of the pairs of B coefficients are of interest, as they correspond to the CP-even and CPodd components of the polarization. The results obtained using the measured coefficients and their covariance TABLE III. Measured coefficients and asymmetries and their total uncertainties. Predicted values from simulation are quoted with a combination of statistical and scale uncertainties, while the NLO calculated values are quoted with their scale uncertainties [3,4]. The NNLO QCD prediction for A jΔϕ ll j , with scale uncertainties, is 0.115 þ0.005 −0.001 [69].   Table VI, and are consistent with the SM predictions. For the coefficients in Table III sensitive to P-and CP-even spin correlations (which are substantial in the SM), we use the NLO calculations to transform the measurements into determinations of f SM , the strength of the given measure of spin correlations relative to the SM prediction. A linear dependence of f SM on the measured coefficient is defined, where f SM ¼ 1 and f SM ¼ 0 correspond to measurements in agreement with the NLO calculations in the presence and absence of spin correlations, respectively. The resulting measurements of f SM are shown in Table VII, where the theoretical scale uncertainty from the transformation is shown as a separate uncertainty. There is a potential correlation between the theoretical scale uncertainty and the scale component of the experimental systematic uncertainty. A similar correlation may exist with the top quark p T systematic uncertainty, owing to its connection to missing higherorder QCD terms. These effects are neglected because their effect on the total uncertainty would be small. The f SM results are also shown in Fig. 13. The results are all consistent with unity, demonstrating the agreement of the measured spin correlation strengths with the SM predictions for all considered combinations of reference axes. FIG. 9. Measured values of the polarization coefficients (circles) and the predictions from POWHEGv2 (triangles), MADGRAPH5_aMC@NLO (inverted triangles), and the NLO calculation [4] (squares). The inner vertical bars on the circles give the statistical uncertainty in the data and the outer bars give the total uncertainty. The numerical measured values with their statistical and systematic uncertainties are given on the right. The vertical bars on the values from simulation represent the combination of statistical and scale uncertainties, while for the calculated values they represent the scale uncertainties. FIG. 10. Measured values of the spin correlation coefficients and asymmetries (circles) and the predictions from POWHEGv2 (triangles), MADGRAPH5_aMC@NLO (inverted triangles), the NLO calculation [3,4] (squares), and the NNLO calculation [69] (cross). The inner vertical bars on the circles give the statistical uncertainty in the data and the outer bars give the total uncertainty. The numerical measured values with their statistical and systematic uncertainties are given on the right. The vertical bars on the values from simulation represent the combination of statistical and scale uncertainties, while for the calculated values they represent the scale uncertainties.
Cross correlation coefficient FIG. 11. Measured values of the cross spin correlation coefficients (circles) and the predictions from POWHEGv2 (triangles), MADGRAPH5_aMC@NLO (inverted triangles), and the NLO calculation [4] (squares). The inner vertical bars on the circles give the statistical uncertainty in the data and the outer bars give the total uncertainty. The numerical measured values with their statistical and systematic uncertainties are given on the right. The vertical bars on the values from simulation represent the combination of statistical and scale uncertainties, while for the calculated values they represent the scale uncertainties.

A. Constraining the top quark CMDM
Analogous to the magnetic dipole moment of an electrically charged particle, the chromomagnetic dipole moment (CMDM) of a color-charged particle in color fields can be defined. In the SM, the intrinsic spin of the top quark and its color charge give it a small CMDM [3,4]. Several BSM models, such as two-Higgs-doublet models (e.g., supersymmetry), technicolor, and top quark compositeness models [71,72], predict an anomalous CMDM, leading to modifications of the tt production rate and spin structure. As a consequence, the measurement of the tt production spin density matrix represents a powerful probe of the top quark CMDM and can be used to search for BSM phenomena.
As in Ref. [12], the effect of an anomalous CMDM on tt production is predicted using an EFT framework in which a fixed set of dimension-six operators is added to the SM Lagrangian [73,74]. The anomalous CMDM of the top quark is a consequence of the O tG operator [71], where y t denotes the Yukawa coupling of the top quark, g S is the strong coupling (g S ¼ 2 ffiffiffiffiffiffiffi ffi πα S p ), Q is the left-handed third-generation quark doublet, σ μν are the Dirac matrices, T a are the Gell-Mann matrices divided by 2, t is the righthanded top quark singlet,φ is the charge-conjugated Higgs doublet field, and G a μν is the gluon field strength tensor. Besides modifying the gtt vertex, O tG also leads to a new ggtt vertex. The contribution due to O tG is parametrized by a dimensionless Wilson coefficient divided by the square of the BSM scale (Λ), assumed to be large compared to the scales typically probed at the LHC. The real part of this Wilson coefficient is denoted as C tG . The imaginary part corresponds to a top quark chromoelectric dipole moment (CEDM), and is assumed to be zero in this section. The top quark CEDM is constrained in Sec. IX B.
To produce predictions for the normalized tt differential cross section, the model of Ref. [71] is implemented in the MADGRAPH5_aMC@NLO generator at NLO in QCD. The setup is similar to that of the MADGRAPH5_aMC@NLO sample introduced in Sec. IV, but without extra partons at the ME level. The RIVET framework [75] is used to apply the object definitions and calculate the spin density matrix observables.
Four observables are chosen to constrain C tG =Λ 2 , corresponding to the four dimensions in Eq. (6), with TABLE V. Summary of the systematic, statistical, and total uncertainties in the extracted tt spin correlation coefficients and asymmetries. An ellipsis (Á Á Á) is shown where the values are <0.0005.

Uncertainty
Source C kk C rr C nn C rk þ C kr C rk − C kr C nr þ C rn C nr − C rn C nk þ C kn C nk − C kn D A lab cos φ A jΔϕ ll j Trigger 0.001 0.001 Á Á Á 0.002 the restriction that they are independent from each other. For example, only two of the observables cos θ k 1 , cos θ r 1 , and cos θ n 1 are independent because they are the direction cosines of the fk;r;ng coordinate system. The C kk , C nn , C rk þ C kr , and D coefficients are all directly sensitive to C tG =Λ 2 [4], and the corresponding four observables (as defined in Table I) are chosen.
A χ 2 minimization technique is used to constrain C tG =Λ 2 , with where data i and pred i ðC tG =Λ 2 Þ are the measured and predicted normalized differential cross sections in the ith of the N bins of the chosen observables, and Cov −1 ij is the (ith, jth) element of the inverse of the data covariance matrix for those N bins. The covariance matrix, corresponding to a subset of the bins illustrated in Fig. 8, accounts for all systematic and statistical uncertainties, as well as the interbin correlations introduced in the unfolding process. In order to break the linear dependencies between the bins of each distribution after normalization to unit area, one bin from each distribution is excluded from the fit, along with the rows and columns associated with it in the covariance matrix. The fit result is independent of the choice of excluded bins.
The χ 2 minimization procedure is performed twice: first including the full contribution from C tG =Λ 2 to the tt cross section, and second including only the contribution that is linear in C tG =Λ 2 , which describes the interference of the O tG amplitudes with those of the SM. In both cases the best-fit value of C tG =Λ 2 is 0.06 TeV −2 , corresponding to a χ 2 =d:o:f: of 8=19. The difference between the two results is negligible, indicating that the value of C tG =Λ 2 is small enough to justify the linear approximation.
Assuming Gaussian probability density functions for the uncertainties in the unfolded data, constraints with confidence levels (C.L.) can be estimated from the values of C tG =Λ 2 for which the Δχ 2 reaches certain values. The Δχ 2 is defined as the change in χ 2 from its minimum value, and is shown as a function of C tG =Λ 2 in Fig. 14. Since the uncertainties in the theoretical predictions do not have a clear frequentist interpretation, they are not included in the confidence intervals. They are estimated separately in Fig. 14 from the maximally positive and negative effects on the best-fit value of C tG =Λ 2 when changing μ R and μ F  13. Measured values of f SM , the strength of the measured spin correlations relative to the SM prediction. The inner vertical bars give the statistical uncertainty, the middle bars give the total experimental uncertainty (statistical and systematic), and the outer bars give the total uncertainty. The numerical measured values with their uncertainties are given on the right.
individually and simultaneously up and down by a factor of 2 in the MADGRAPH5_aMC@NLO predictions.
The resulting constraint at the 95% C.L. is −0.10 < C tG =Λ 2 < 0.22 TeV −2 . In Ref. [71], a 95% C.L. constraint of −0.42 < C tG =Λ 2 < 0.30 TeV −2 was derived using NLO predictions for the contributions from O tG to the total tt cross section, combined with ATLAS and CMS measurements at ffiffi ffi s p ¼ 8 TeV, as well as −0.32 < C tG =Λ 2 < 0.73 TeV −2 , using Fermilab Tevatron results. From a measurement of the absolute tt differential cross section as a function of jΔϕ ll j at the particle level, CMS determined −0.06 < C tG =Λ 2 < 0.41 TeV −2 at the 95% C.L. [12]. The results presented here are consistent with and improve on these previous limits. Compared to Ref.
[12], the sensitivity to C tG =Λ 2 in this analysis is improved by 30% and the theoretical uncertainties are substantially smaller.

B. Constraining anomalous couplings
The top quark anomalous CMDM operator is just one of the 11 independent dimension-six operators relevant for hadronic tt production [4]. The normalized differential cross sections measured in Sec. VIII A are sensitive to ten of these operators, and each can be constrained using a fit similar to that in Sec. IX A. However, in the absence of a consistent simulation of all these operators compatible with the NLO QCD predictions in the MADGRAPH5_aMC@NLO generator [71], we instead use the known functional forms of Eqs. (8)-(10) to fit the data. We use the calculations from Ref. [4] to determine the coefficients and their dependence on the contributions from the different operators. The calculations for the NLO SM part are the same as those introduced in Sec. VIII A. For the contributions from the operators, only tree-level interference terms with the QCD amplitudes in the linear approximation are considered [4].
The anomalous couplings associated with the 11 operators are listed in Table VIII, with a brief description of their properties. Unlike the Wilson coefficient C tG considered in Sec. IX A, the anomalous couplings apply to operators in their form after spontaneous symmetry breaking [4]. The couplingsμ t andd t represent the top quark anomalous CMDM and CEDM, respectively, and there are two further CP-odd operators involving two top quarks and up to three gluons (with couplingsĉ −− andĉ −þ ). The operators associated with the remaining couplings represent CP-even four-quark interactions, with weak isospin quantum numbers either 0 or 1. The operators are described in detail in Ref. [4].
The normalized differential cross sections measured in Sec. VIII A are sensitive to all the anomalous couplings given in Table VIII exceptĉ AA , which is constrained by measurements of the tt charge asymmetry [4]. Using the same fitting procedure as in Sec. IX A, we set a limit on each coupling, setting the other couplings to zero. The 95% C.L. limits are given in Table IX, and the measured values and uncertainties are listed and displayed in Fig. 15. Theoretical uncertainties are estimated from the simultaneous variation of μ R and μ F up and down by a factor of 2. Limits are given for the combination of couplingŝ c 1 −ĉ 2 þĉ 3 rather thanĉ 2 alone because this is the combination of couplings to which the measurements are The Δχ 2 values from the fit to the data as a function of C tG =Λ 2 . The solid line is the result of the nominal fit, and the dotted and dashed lines show the most-positive and mostnegative shifts in the best-fit C tG =Λ 2 , respectively, when the theoretical inputs are allowed to vary within their uncertainties. The vertical line denotes the best-fit value from the nominal fit, and the inner and outer areas indicate the 68 and 95% C.L., respectively. TABLE VIII. Anomalous couplings associated with the dimension-six operators relevant for hadronic tt production, the operator type of the effective interaction vertex they represent, and their P and CP symmetry properties. It is not possible to combine the isospin-1 operators such that they have definite properties with respect to C and P [4].

Coupling
Operator type Symmetry propertieŝ μ t 2 quarks plus gluon(s) P-even, CP-even d t 2 quarks plus gluon(s) P-odd, CP-odd c −− 2 quarks plus gluon(s) P-odd, CP-odd c −þ 2 quarks plus gluon(s) P-even, CP-odd c VV 4 quarks (weak isospin 0) P-even, CP-even c VA 4 quarks (weak isospin 0) P-odd, CP-even c AV 4 quarks (weak isospin 0) P-odd, CP-even c AA 4 quarks (weak isospin 0) P-even, CP-even c 1 4 quarks (weak isospin 1) CP-even c 2 4 quarks (weak isospin 1) CP-even c 3 4 quarks (weak isospin 1) CP-even directly sensitive [4]. The strongest constraints are found for the operators probed in the gg initial state. The fourquark operators with isospin 0 are more constrained than those with isospin 1, where contributions from the up and down quark qq initial states have opposite signs and similar magnitudes [4,16]. We also consider the simultaneous fitting of multiple couplings. We find that the pairs of four-quark couplings ðĉ VV ;ĉ 1 Þ, ðĉ VA ;ĉ 3 Þ, and ðĉ AV ;ĉ 1 −ĉ 2 þĉ 3 Þ cannot be simultaneously constrained because their predicted effects on the measured distributions can approximately cancel each other. The constraints on the other couplings are found to be independent, and therefore sufficiently characterized by the results of Table IX, with the exception of three combinations of couplings for which we derive two-dimensional 68 and 95% C.L. limits, shown in Fig. 16.
For a direct comparison with the top quark CMDM results of Sec. IX A, we use the relationship C tG =Λ 2 ¼μ t =ð2m t 2 Þ. Taking the result forμ t from Table IX, we find a central value of C tG =Λ 2 ¼ −0.09 TeV −2 , with −0.24 < C tG =Λ 2 < 0.07 TeV −2 at the 95% C.L. The sensitivity to C tG =Λ 2 (determined from the width of the confidence interval) is the same as that found in Sec. IX A, which suggests that the tree-level calculation of the interference terms, using the linear approximation, is adequate for C tG =Λ 2 . The difference in central value is attributable to the difference in the SM predictions for the coefficients in the NLO calculations and the MADGRAPH5_aMC@NLO simulation. Since the SM prediction is of greater accuracy in the NLO calculations (which include EW corrections), we quote the C tG =Λ 2 result of this section as the nominal result of the analysis.
In a similar way,d t is related to the imaginary part of the Wilson coefficient of the O tG operator C I tG , and we find a constraint at the 95% C.L. of −0.33<C I tG =Λ 2 <0.20TeV −2 , with a central value of C I tG =Λ 2 ¼ −0.07 TeV −2 . This represents a substantial improvement over existing direct constraints on the top quark CEDM [76,77], but it is still significantly weaker than the indirect constraint of jC I tG =Λ 2 j < 0.007 TeV −2 [78] derived from the experimental limit on the neutron electric dipole moment [79,80].
Analogous to the magnetic and electric dipole moments, μ t andd t can be expressed in terms of the dimensionful parameters C 5 and D 5 , which are related to the former by a factor of 1=m t [81]. In this parametrization, we find constraints at the 95% C.L. of ð−1.6 <C 5 < 0.5Þ× 10 −18 g S cm and ð−2.3 < D 5 < 1.4Þ × 10 −18 g S cm.  Table VIII, derived by fitting the distributions measured in Sec. VIII A and setting the other anomalous couplings to zero. The confidence intervals include only the experimental uncertainties as in Sec. IX A. The theoretical uncertainties, the χ 2 values (d:o:f: ¼ 19), and the distributions used in each fit are given in the last three columns. For conciseness, the distributions are labeled by their associated coefficients (as defined in Table I). An ellipsis (Á Á Á) is shown where the uncertainties are <0.0005. 95% C.L.

X. SUMMARY
Measurements of the top quark polarization and tt spin correlations have been presented, probing all of the independent coefficients of the top quark spin-dependent parts of the tt production density matrix for the first time in proton-proton collisions at ffiffi ffi s p ¼ 13 TeV. Each coefficient was extracted from a normalized differential cross section, unfolded to the parton level and extrapolated to the full phase space. The measurements were made using a data sample of events containing two oppositely charged leptons (e þ e − , e AE μ ∓ , or μ þ μ − ) and two or more jets, of which at least one was identified as coming from the hadronization of a bottom quark. The data were recorded by the CMS experiment in 2016 and correspond to an integrated luminosity of 35.9 fb −1 . The measured normalized differential cross sections and coefficients were compared with standard model predictions from simulations with NLO accuracy in QCD and from NLO QCD calculations including electroweak corrections. The measured distribution of jΔϕ ll j, the absolute value of the difference in azimuthal angle between the two "Theory unc. up" refers to the fit value when μ R and μ F are simultaneously increased by a factor of 2, and "theory unc. down" when they are decreased by the same factor. leptons in the laboratory frame, was additionally compared with a NNLO QCD prediction. All of the measurements were found to be consistent with the expectations of the standard model. The distribution of cos φ, equivalent to the dot product of the two lepton directions measured in their parent top quark and antiquark rest frames, is most sensitive to the presence of spin correlations, with a relative uncertainty below 5%. Statistical and systematic covariance matrices were provided for the set of all measured bins, and were used in simultaneous fits to constrain the contributions from ten dimension-six effective operators. Two of these operators represent the anomalous chromomagnetic and chromoelectric dipole moments of the top quark, and constraints on their Wilson coefficients of −0.24 <C tG =Λ 2 < 0.07 TeV −2 and −0.33 < C I tG =Λ 2 < 0.20 TeV −2 , respectively, were obtained at the 95% confidence level. This constitutes a substantial improvement over previous direct constraints.

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 centers 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