Measurement of the top quark Yukawa coupling from $\mathrm{t\bar{t}}$ kinematic distributions in the dilepton final state in proton-proton collisions at $\sqrt{s}=$ 13 TeV

A measurement of the Higgs boson Yukawa coupling to the top quark is presented using proton-proton collision data at $\sqrt{s}=$ 13 TeV, corresponding to an integrated luminosity of 137 fb$^{-1}$, recorded with the CMS detector. The coupling strength with respect to the standard model value, $Y_\mathrm{t}$, is determined from kinematic distributions in $\mathrm{t\bar{t}}$ final states containing ee, $\mu\mu$, or e$\mu$ pairs. Variations of the Yukawa coupling strength lead to modified distributions for $\mathrm{t\bar{t}}$ production. In particular, the distributions of the mass of the $\mathrm{t\bar{t}}$ system and the rapidity difference of the top quark and antiquark are sensitive to the value of $Y_\mathrm{t}$. The measurement yields a best fit value of $Y_\mathrm{t}=$ 1.16$^{+0.24}_{-0.35}$, bounding $Y_\mathrm{t}\lt$ 1.54 at a 95% confidence level.


Introduction
Since the discovery of the Higgs boson in 2012 [1,2], one of the main goals of the CERN LHC program has been to study in detail the properties of this new particle. In the standard model (SM), all fermions acquire their mass through the interaction with the Higgs field. More specifically, the mass of a given fermion, m f , arises from a Yukawa interaction with coupling strength where v is the vacuum expectation value of the Higgs field. Among all such couplings, the top quark Yukawa coupling is of particular interest. It is not only the largest, but also remarkably close to unity. Given the measured top quark mass [3,4], the mass-Yukawa coupling relation implies a value of the Yukawa coupling g SM t ≈ 0.99 when evaluated near the energy scale of m t . Physics beyond the SM, such as two Higgs doublet and composite Higgs boson models, introduce modified couplings that alter the interaction between the top quark and the Higgs field [5,6]. This makes the interaction of the Higgs boson with the top quark one of the most interesting features of the Higgs field to study at the LHC today, especially because it is experimentally accessible through multiple avenues, both direct and indirect.
For the purpose of this measurement, we define for the top quark the parameter Y t = g t /g SM t , which is equivalent to the modifier κ t introduced in the κ-framework [7]. We consider only the case where Y t ≥ 0, though certain specific techniques are sensitive also to the sign of the Yukawa coupling (for example, Ref. [8]). Recent efforts have had notable success in directly probing g t via the production of a Higgs boson in association with a top quark pair (ttH) [9, 10]. Currently, the most precise determination comes from the κ-framework fit in Ref. [11], which yields Y t = 0.98 ± 0.14 by combining information from several Higgs boson production and decay channels. These measurements, however, fold in assumptions of the SM branching fractions via Higgs couplings to other particles. Another way to constrain g t , which does not depend on these couplings, was presented in the search for four top quark production in Ref. [12], yielding a limit of Y t < 1.7 at a 95% confidence level (CL). However, it is also possible to constrain g t indirectly using the kinematic distributions of reconstructed tt pair events, a technique that has been recently used by CMS to derive a similar limit of Y t < 1.67 at 95% CL in the lepton+jets tt decay channel [13]. The measurement presented in this paper follows this last approach, but in the dilepton final state.
Current commonly used Monte Carlo (MC) simulations of tt production include next-to-leadingorder (NLO) precision in perturbative quantum chromodynamics (QCD). Subleading-order corrections arise from including electroweak (EW) terms in the perturbative expansion of the strong coupling α S and the EW coupling α. Such terms begin to noticeably affect the cross section only at loop-induced order α 2 S α, and are typically not included in current MC simulation. While these terms have a very small effect on the total cross section, they can alter the shape of kinematic distributions to a measurable extent. Such changes become more noticeable if the Yukawa coupling affecting the loop correction ( Fig. 1) is anomalously large. Therefore, these corrections are of particular interest in deriving upper limits on g t . For example, the distribution of the invariant mass of the tt system, M tt , will be affected significantly by varying Y t . Doubling the value of Y t can alter the M tt distribution by about 9% near the tt production threshold, as described in Ref. [14]. Another variable sensitive to the value of Y t is the difference in rapidity between the top quark and antiquark, ∆y tt = y(t) − y(t). In tt production, M tt and ∆y tt are proxies for the Mandelstam kinematic variables s and t, respectively, which span the event phase space and can thus be used to include the EW corrections in previously generated event samples via reweighting. The effects of these corrections are shown for differential cross sections of M tt and ∆y tt in Fig. 2. These are computed by reweighting simulated tt events at generator level using predictions from the HATHOR software package [15].  After calculating the dependence of these corrections on Y t , a measurement is performed. We use events in the dilepton final state (ee, µµ, or eµ), for which this type of measurement has not yet been performed. While this decay channel has a smaller branching fraction than the lepton+jets channel studied in Ref.
[13], it has lower backgrounds due to the presence of two final-state high-p T leptons. However, two neutrinos are also expected in this final state, which escape detection and pose challenges in the kinematic reconstruction. For this reason, we do not perform a full kinematic reconstruction as was done in the previous measurement in the lepton+jets channel. This measurement also utilizes a much larger data set with an integrated luminosity of 137 fb −1 collected during Run 2 at the LHC from 2016 to 2018, allowing us to achieve comparable precision to that in Ref.
[13] for a decay channel with a much lower branching fraction.

Simulation of top quark pair production and backgrounds
The production of tt events is simulated at the matrix-element (ME) level with NLO QCD precision, using the POWHEG 2.0 generator [19][20][21][22]. The calculation is performed with the renormalization and factorization scales, µ R and µ F , set to the transverse top quark mass, m T = √ m 2 t + p 2 T , where p T is the transverse momentum of the top quark and the quantity is evaluated in the tt rest frame. The default value of m t is set to 172.5 GeV. The ME calculations obtained from POWHEG are combined with the parton shower simulation from PYTHIA 8.219 [23], using the underlying-event tune M2T4 [24] to simulate data taken in 2016, and PYTHIA 8.226 using the tune CP5 [25] to simulate data taken in 2017 and 2018. The parton distribution function (PDF) set NNPDF3.0 at NLO [26] is used for 2016 and updated to NNPDF3.1 [27] at next-to-NLO (NNLO) for 2017 and 2018. These samples are normalized to a tt cross section calculated at NNLO in QCD including resummation of next-to-next-to-leading logarithmic (NNLL) soft gluon terms with TOP++ 2.0 [28]. The calculation uses the PDF4LHC prescription [29] with the MSTW2008 NNLO [30,31], CT10 NNLO [32,33] and NNPDF2.3 [34] PDF sets used to generate an envelope of uncertainty with the midpoint of the envelope used for the central predictions. The PDF uncertainty is then summed in quadrature with the scale uncertainty to arrive at an overall uncertainty of ≈5% on the nominal value of 832 pb. The shape effects associated with the PDF uncertainty are considered separately in Section 7.
A high purity of tt events can be obtained in the dilepton channel, as shown in Section 4. A small contamination is expected to result from background processes, which are modeled by simulation. In particular, we account for dilepton production due to Drell-Yan type processes and single top quark production. Other SM processes, such as W boson production, were investigated and found to have negligible contributions. Diboson production is also included, although its expected contribution is minute due to the small total cross section of the process.
About 1% of the events identified as tt dilepton decays are misidentified tt lepton+jets decays. EW corrections are applied to all tt events, even misidentified ones, so their kinematic distributions remain dependent on Y t . Thus, these events are still considered as signal, even though their contribution to the measurement sensitivity is greatly diminished relative to dilepton events.
Single top quark events are simulated at NLO with POWHEG in combination with PYTHIA, while diboson events are simulated with PYTHIA at leading-order (LO) QCD precision. Drell-Yan production is simulated at LO using MADGRAPH5 aMC@NLO version 2.2.2 for 2016 and version 2.2.4 for 2017 onwards [35], with up to 4 additional partons, interfaced to PYTHIA using the MLM matching algorithm [36,37].
The detector response to all simulated events is modeled with the GEANT4 software toolkit [38]. In addition, the effects of multiple proton-proton interactions per event are included in simulations and the distribution of these pileup interactions is reweighted to the vertex multiplicity distribution in the data.

Simulation of electroweak corrections
Contributions to the top quark pair production arising from QCD+EW diagrams are evaluated using the HATHOR package [15], which is used to compute a double-differential cross section as a function of M tt and ∆y tt including LO QCD diagrams and certain EW diagrams of order α 2 S α. These diagrams involve massive boson exchange and examples are shown in Fig. 1. The contributions from photon-mediated interactions are not included. Contributions from diagrams involving virtual photon exchange should not be assessed individually, as they are partially cancelled not only by real emission diagrams but also by contributions from gγ production [39]. A complete assessment would require the modeling of photon content within protons. This was not performed here, as the net effect is fairly small. For example, Ref. [39] cites a 1% effect from photon-mediated contributions to the tt cross section at the LHC with detector-based kinematic cuts. Thus, we include only diagrams involving massive vector and scalar boson interactions, which are the dominant EW diagrams at this order.
The ratio of this double-differential cross section is evaluated with respect to the LO QCD computation, in order to obtain a multiplicative weight correction w(M tt , ∆y tt ). Applying this weight at parton level to MC samples produced at NLO QCD approximates the inclusion of EW corrections in the simulation. This multiplicative approach to including EW corrections was used previously in Ref. [13], and has the benefit of approximating the inclusion of diagrams at order O(α 3 S α). Because EW corrections factorize in some kinematic regimes, this is a better-motivated approach than the alternative additive approach, in which one adds the fixed-order result at order O(α 2 S α) while ignoring all potential contributions of order O(α 3 S α) (see Ref. [40] for a more detailed discussion). In other words, the additive approach applies the EW correction factor only to the proportion of POWHEG events present at LO QCD, while any interplay between EW corrections and higher-order QCD simulation is ignored. Although the multiplicative approach is clearly favored, neither approach can account for the effects of two-loop contributions near the tt production threshold. To account for the lack of knowledge of such terms, we take the difference between the two predictions as a modeling uncertainty in this regime, as suggested in Ref. [14]. The estimation of this uncertainty is discussed further in Section 6.
The EW correction weights are calculated for discrete integer values of Y t = 0, 1, . . . , 5. Since the dependence of the production rate on Y t is exactly quadratic, these discrete values are sufficient to parametrize event yields as a continuous function of Y t (as discussed in Section 6). This allows us to measure which value of Y t best describes the data.

Event and object selection
Events are selected using single-electron or single-muon triggers. Data taking at the LHC was interrupted by technical stops at the end of each year, leading to some changes in configuration and modeling between 2016, 2017, and 2018. For events selected by the single-electron trigger, we require a trigger p T threshold of 27 GeV with the exception of 2018, where a threshold of 32 GeV is used. In the case of the single-muon trigger, we select events with a trigger p T threshold of 24 GeV, which is raised to 27 GeV only for 2017 due to high event rates.
We ensure that all electrons and muons are within the silicon tracker coverage by requiring a pseudorapidity |η| < 2.4. To operate well above the trigger threshold, we then require at least one isolated electron or muon reconstructed with p T > 30 GeV, except in 2018, where we require leading p T electrons to have p T > 34 GeV in accordance with the trigger threshold. The same lepton isolation criteria described in Ref. [41] are used. After selecting the leading p T lepton, a second isolated electron or muon with p T > 20 GeV is required. Events with three or more isolated leptons with p T > 15 GeV are discarded.
Jets are clustered from PF objects via the anti-k T algorithm [42,43] with a distance parameter of 0.4. The jet momentum is calculated as the vectorial sum of the momenta of its constituents. Corrections to the jet energy are derived as a function of jet p T and η in simulation and improved by measurements of energy balance in data [44]. We select jets with |η| < 2.4 and p T > 30 GeV.
Jets originating from b quarks are identified using the DeepCSV algorithm [45]. The algorithm provides three working points: loose, medium, and tight, in order of decreasing efficiency and increasing purity. The b identification efficiencies (and light quark misidentification rates) are 84 (11)%, 68 (1.1)%, and 50 (0.1)%, respectively. For an initial selection, we consider events with a minimum of two b jet candidates passing the loose working point of the algorithm. When applied to simulated tt dilepton decays, we find that this initial selection of b jets will correctly include both b jets originating from top quark decays in 87% of events. In around 9% of simulated tt dilepton events passing this initial selection, there are more than two jets passing the loose working point, leading to an ambiguity in jet assignment. If such events have exactly two jets passing a higher working point (medium or tight), then those two jets are considered the viable candidates for b jets originating from a top quark decay, and the ambiguity is resolved without using kinematic properties of the event. The small fraction of events with more than two viable b jet candidates, making up 4% of the initially selected tt dilepton events, are discarded. After this selection procedure, each event remaining in the sample has exactly two b jet candidates, which together are correctly identified in 85% of simulated tt dilepton events.
In order to remove Drell-Yan background events in the ee and µµ channels, we reject events in which the two leptons have an invariant mass below 50 GeV or within 10 GeV around the Z boson mass of 91.2 GeV.
The missing transverse momentum vector ( p miss T ), defined as the negative vector sum of all transverse momenta, is generally of large magnitude in dilepton decays because of the two undetected final-state neutrinos. To further aid in removing Drell-Yan events, we impose an additional selection requirement on the magnitude of the missing transverse momentum, requiring p miss T > 30 GeV in all events with ee or µµ in the final state.
The breakdown of expected signal and background yields, summed over the three channels (ee, µµ, eµ), is shown by year in Table 1. The Drell-Yan background is estimated to be about 2%. Single top quark production accounts for roughly another 2% of the estimated sample composition. Table 1: Simulated signal, background, and data event yields for each of the three years and their combination. The rightmost column shows the fraction of each component relative to the total simulated sample yield across the full data set. The statistical uncertainty in the simulated event counts is given.

Event reconstruction
The EW corrections are calculated based on M tt and ∆y tt . However, to evaluate these quantities it is necessary to reconstruct the full kinematic properties of the tt system, including the two undetected neutrinos. While it is possible to completely reconstruct the neutrino momenta in the on-shell approximation, such a reconstruction is highly sensitive to p miss T , which introduces large resolution effects and additional systematic uncertainties. We observe that using the proxy variables M bb = M(b + b + + ) and |∆y b b | = |y(b + ) − y(b + )|, where represents a final-state electron or muon, results in a more precise measurement.
Unlike M bb , the accurate reconstruction of |∆y b b | requires that each of the two b jets is matched to the correct lepton, i.e., both originating from the same top quark decay. In order to make this pairing, we utilize the information from the kinematic constraints governing the neutrino momenta.
If one assumes the top quarks and W bosons to be on-shell, the neutrino momenta are constrained by a set of quadratic equations arising from the conservation of four-momentum at each vertex. We refer to these kinematic equations, collectively, as the mass constraint. The mass constraint for each top quark decay results in a continuum of possible solutions for neutrino momenta, which geometrically can be presented as an intersection of ellipsoids in threedimensional momentum-space [46]. For certain values of input momenta of b jets and leptons these ellipsoids do not intersect at all, such that the quadratic equations have no real solution. In these scenarios, the mass constraint cannot be satisfied.
In cases where the mass constraint can be satisfied, one could also constrain p miss T in the event to equal the p T sum of the two undetected neutrinos. We call this the p miss T constraint. This constraint reduces the remaining solutions to a discrete set, containing either two or four possibilities that fully specify the momenta of both neutrinos. Similar to the case of the mass constraint, there are some values of the input parameters for which the p miss T constraint cannot be satisfied.
When looking at simulated events where both b jets are correctly reconstructed and paired, we find that the mass constraint can be satisfied in 96% of all cases, while the mass and p miss T constraints can be simultaneously satisfied in 55% of cases. In contrast, if the b jets are correctly reconstructed but incorrectly paired to leptons, the mass constraint can be satisfied in only 23% of cases, while both mass and p miss T constraints can be met in only 18% of cases.
Pairings with no solution to the mass constraint are thus frequently incorrect. When the mass constraint can be satisfied, pairings with a solution to the p miss T constraint are more likely to be correct. This information is used as part of the pairing procedure, which has three steps.
1. The mass constraint is checked for both possible pairings. If only one pairing is found to satisfy the mass constraint, that pairing is used. If both pairings fail to satisfy the mass constraint, the event is discarded. If both pairings satisfy the mass constraint, we check the p miss T constraint.
2. If only one pairing allows for the p miss T constraint while the other does not, the pairing yielding an exact solution to the p miss T constraint is used.
3. If the kinematic variables of the neutrinos do not suggest a clear pairing, the b jets, b 1 and b 2 , are paired with the leptons ( , ) by minimizing the quantity and ϕ is the azimuthal angle in the transverse plane.
In simulation, this procedure discards 7% of the signal sample, targeting events which generally involve an improperly assigned or misidentified b jet (at a rate of 72%). This raises the fraction of events that successfully identify both b jets from a top quark decay to 89% in simulation. After these steps, we obtain the correct b jet pairing in 82% of simulated dilepton tt events for which both b jets originating from top quark decays were correctly identified, and thus 73% of simulated dilepton tt events overall.
The sensitivity of our chosen kinematic variables to Y t , before and after reconstruction, is shown in Fig. 3. We see that, in the chosen proxy variables, not much sensitivity is lost in the reconstruction process. This is especially true for the proxy mass observable, M bb , providing an advantage over M tt , which cannot be reconstructed as accurately.

Comparison between data and simulation
Comparisons between data and simulation are shown in Fig. 4, where tt events are broken into four categories: events with correctly identified leptons and jets in which jets are correctly assigned (tt correct jets), events with correctly identified leptons and jets in which jets are incorrectly assigned (tt swapped jets), events with correctly identified leptons where the two b jets originating from top quark decays are not identified correctly (tt wrong jets), and lastly events where the identified leptons are not those from W decay vertices (tt wrong leptons). The majority of events in the last category are tt dilepton decays where a W boson decay produces a τ lepton which itself decays leptonically, with a small fraction being misidentified decays in the lepton+jets channel (1% of the total tt signal). Though all tt events are subject to EW corrections and thus considered as signal, the sensitivity of the reconstructed kinematic variables is generally decreasing among the four categories.
Various observations can be made from Fig. 4. The agreement between data and simulation appears generally to be within the total uncertainty (discussed further in Section 7), and the     Figure 3: The ratio of kinematic distributions with EW corrections (evaluated for various values of Y t ) to the SM kinematic distribution (Y t =1) is shown, demonstrating the sensitivity of these distributions to the Yukawa coupling. The plots on the left show the information at the generator level, while the plots on the right are obtained from reconstructed events. The axis scale is kept the same for the sake of comparison.
small overall background rate is apparent. Most events are seen to be associated with zero or one additional jet (beyond the two b jets). The effect of the p miss T selection requirement can be seen, removing events in the ee and µµ final states in a regime with high Drell-Yan background rates. Single top quark production background rates are seen to vary less steeply as a function of p miss T . Looking at the leading lepton p T , we see that the additional use of a dilepton trigger would not yield a substantial increase in sensitivity.
A slope is apparent in the ratio of data to the MC prediction in the p T distributions of leptons and b jets. The trends may be related to a previously observed feature of the nominal POWHEG+PYTHIA simulation, in which a harder top quark p T distribution is observed than in data (as discussed, e.g., in Ref. [41]). This behavior is the subject of much discussion in the top quark physics community, so we remark on it in this paper despite the fact that we are primarily concerned with other kinematic variables. Fixed-order NNLO calculations are available that generally show a softer top quark p T spectrum than in the POWHEG+PYTHIA simulation, which could be seen as evidence that the discrepancy arises from mismodeling in simulation. However, the modeling of M tt does not appear to suffer such issues [41], and we see no evidence that the kinematic variables used in this measurement are not well-described within the included modeling uncertainties. Further discussion of the top quark p T spectrum in POWHEG+PYTHIA and its relation to fixed-order NNLO calculations can be found in Section 7.   The uncertainty bands are derived by varying each uncertainty source up and down by one standard deviation (as described in Section 7) and summing the effects in quadrature. The signal simulation is divided into the following categories: events with correctly identified leptons and jets in which jets are correctly assigned (tt correct jets), events with correctly identified leptons and jets in which jets are incorrectly assigned (tt swapped jets), events with correctly identified leptons where the two b jets originating from top decays are not identified correctly (tt wrong jets), and lastly events where the identified leptons are not those from W boson decay vertices (tt wrong leptons). The lower panels show the ratio of data to the simulated events in each bin, with total uncertainty bands drawn around the nominal expected bin content.

Measurement strategy and statistical methods
After reconstruction, events are binned coarsely in |∆y b b | and more finely in M bb . The binning is chosen to ensure each bin in each data-taking year contains at least on the order of 10 000 events, as seen in Fig. 5, leading to a low statistical uncertainty and improved uncertainty estimation.
In each bin, the expected yield is parametrized as a function of Y t . The effect is exactly quadratic, as a consequence of the order at which EW corrections are evaluated. We perform a quadratic fit to extrapolate the effect of the EW corrections on a given bin as a continuous function of Y t (Fig. 6). This correction for each bin can be applied as a rate parameter R EW affecting the expected bin content.   Figure 5: The pre-fit agreement between data and MC simulation in the final kinematic binning. The solid lines divide the three data-taking periods, while the dashed lines divide the two |∆y b b | bins in each data-taking period, with M bb bin ranges displayed on the x axis. The lower panel shows the ratio of data to the simulated events in each bin, with total uncertainty bands drawn around the nominal expected bin content, obtained by summing the contributions of all uncertainty sources in quadrature.
We construct a likelihood function L, where φ and {θ i } are the suite of nuisance parameters associated with individual sources of systematic uncertainty. The distributions p(φ) and p(θ i ) are penalty terms which assign probability distributions that encode the likelihood the parameters vary from their prior values, as discussed further below. Each bin has an individual Poisson likelihood distribution, describing the probability of a bin content to vary from statistical fluctuations. Here n bin obs is the total observed bin count, with the expected bin count being the sum of the predicted signal yield s bin and background yield b bin . The number of expected signal events is modified by the additional rate parameter R EW , which depends on the Yukawa coupling ratio Y t and a special nuisance parameter φ that encodes the uncertainty associated with the multiplicative application of EW corrections derived at order O(α 2 S α). The full expression for the rate R bin EW , including this uncertainty term in the bins near the tt production threshold, is given by  (4) In the nominal case, we have R bin EW (Y t ) = 1 + δ EW (Y t ). Intuitively, δ EW represents the marginal effect of EW corrections included in HATHOR relative to the LO QCD calculation, while δ QCD represents the marginal effect of higher-order terms included in the POWHEG sample relative to the LO QCD calculation. The multiplicative approach to including EW corrections assumes that these two corrections factorize. The quantity δ bin QCD δ bin EW represents the cross term arising from the difference in multiplicative and additive approaches. The Gaussian distributed nuisance parameter φ modulates the uncertainty generated by this cross term. We note that the uncertainty in the EW corrections is unique because it depends on the value of Y t at which the EW corrections are evaluated. Thus, it is described by its own term and nuisance parameter φ, separate from other systematic uncertainties. For bins away from the threshold where EW corrections decrease as a function of Y t , we do not include this uncertainty. These bins do not contribute much sensitivity to the measurement, and enter a kinematic regime in which this method of uncertainty estimation is no longer meaningful. At the large values of the Mandelstam variable s that correspond to these bins, the dominant terms contributing to δ EW are Sudakov logarithms resulting from W and Z boson exchange. These terms factorize well and do not contribute to the uncertainty we wish to model [40].
Each nuisance parameter θ j corresponding to an overall normalization uncertainty, such as the uncertainty in the integrated luminosity or in cross section values, is assumed to follow a lognormal distribution p(θ j ). Uncertainties with shape effects associated to nuisance parameters {θ i } are handled by generating up and down variations of the bin content s bin for each θ i . These variations result from changing the underlying theoretical/experimental sources, which are outlined in Section 7, usually by one standard deviation (σ) based on the uncertainty in our best estimates. These up and down variations are then enforced to correspond to the bin modifiers associated with θ i = ±1, while θ i = 0 corresponds to the nominal estimate. The nuisance parameter θ i is then taken to follow a Gaussian distribution p(θ i ) with mean µ = 0 and variance σ 2 = 1 in the likelihood. The collection of bin modifiers for these up and down variations are referred to as templates, with examples shown in Section 8. A vertical template morphing is applied to alter the shape as a function of the underlying nuisance parameter θ i , where in each bin the modifier is interpolated as a sixth-order polynomial spline for values of θ i ∈ [−1, 1] and linearly outside of that region, assuring that s bin (θ i ) remains continuous and twice differentiable.
The measurement of Y t is then performed via a profile likelihood scan, as described in Ref. [47]. By repeating a maximum likelihood fit over a fine array of fixed values of Y t and comparing to the likelihood at the best fit value, we can use the properties of the maximum likelihood test statistic to evaluate intervals at 68% and 95% CL around the best fit value.

Experimental and theoretical uncertainties 7.1 Sources of uncertainty
The list of uncertainties considered is very similar to that of the previous measurement presented in Ref. [13]. The main differences are the lack of QCD multijet background and the use of data from the full Run 2 data-taking period. Full or partial correlations are imposed on the underlying uncertainty sources between data-taking periods where appropriate, as discussed further in Section 7.2. Uncertainties that do not alter the shape of the final distribution are treated as normalization uncertainties, while all others are treated as shape uncertainties on the binned data. Shape effects are considered for the distributions of tt events only, as the contribution of background events is small. Correlations of the uncertainties between different data-taking periods are treated on a case-by-case basis. Because the measurement is more sensitive to shape effects than normalization effects, the uncertainties with the largest magnitude do not necessarily have the largest impact on the measurement sensitivity. By repeating the measurement with nuisance parameters frozen at their post-fit values, we evaluate what fraction of the measured uncertainty on Y t is associated to certain nuisance parameters.
The dominant experimental uncertainty in this analysis comes from the calibration of the detector jet energy response. Corrections to the reconstructed jet energies are applied as a function of p T and η. We follow the standard approach outlined in Ref. [44] to consider 26 separate uncertainties that are typically involved in determining these calibrations. In this approach, the uncertainty in the resolution of the jet reconstruction is also considered in addition to the energy response. The effect of these uncertainties is propagated to the reconstruction of p miss T . These effects account for approximately 7% of the total uncertainty on Y t in the final measurement.
Other experimental sources of uncertainty are comparatively minor. The overall uncertainty in the integrated luminosity of 2.5, 2.3, and 2.5% is included as a normalization uncertainty applied to all signal and background events in 2016, 2017, and 2018, respectively [48-50]. The uncertainty in the number of pileup events included in simulation is assessed by varying the inelastic cross section, 69.2 mb, by 4.6% [51].
Efficiencies in b jet identification and misidentification are corrected to match data [45]. While this source is treated as a shape effect, the uncertainty manifests approximately as an overall normalization effect on the signal of around 3%, and contributes only about 1% of the final uncertainty on Y t .
Similarly, scale factors are applied in bins of p T and η to correct simulated efficiencies of lepton reconstruction, identification, isolation, and triggers to match data. These are derived from a fit using the tag-and-probe method using Z boson decays [52][53][54]. This fit accounts for the uncertainty from the limited number of events in the data sample as well as differences in performance based on the jet multiplicity. Overall, the effect is assessed to be below 2%.
As a standard technique to estimate the contributions of higher-order QCD terms at the ME level, the renormalization scale µ R and factorization scale µ F are each varied up and down in the POWHEG simulation by a factor of 2. Templates are generated for the individual variation of µ R and µ F , as well as an additional template for the simultaneous variation of the two scales together (up and down), leading to three separate shape uncertainties in total. Since an NNLO tt cross section is already used to improve the normalization of the MC simulation, the normalization effect induced by the scale variations is overestimated. As we include a separate uncertainty on the cross section normalization, the overall normalization effect is therefore removed entirely from the scale variation templates, which are normalized to the nominal sample. The resulting shape effects remain significant and these are among the limiting uncertainties in the fit, contributing about 7% of the total measurement uncertainty.
A 5% normalization uncertainty is assumed in the tt cross section, which covers expected contributions from the higher-order terms not included in the NNLO+NNLL cross section calculation [28], giving a more realistic normalization uncertainty than the variation of µ R and µ F in POWHEG. The backgrounds in this analysis are small enough (≈2% sample composition each) that we do not generate templates for their response to individual systematic uncertainties. A 15% normalization uncertainty is included on single top quark MC samples, which covers the expected ME scale variation and the jet energy correction uncertainties associated with these samples. The Drell-Yan and diboson MC samples are assigned a 30% normalization uncertainty, to cover the larger ME scale variation uncertainties associated with these LO simulations. The background normalizations can alter slightly the expected shape of the data, but are not among the most impactful uncertainties.
We include an uncertainty in the EW corrections, based on our methods for generating and applying these additional terms, as outlined in Sections 3 and 6. Like the scale variations, this uncertainty is designed to cover higher-order effects at the ME level, specifically those arising from diagrams of order α 3 S α. It places an uncertainty on R bin EW of 10-40% in the applicable bins, which translates to a small overall uncertainty in bin rate unless the corrections are evaluated at a value of Y t far from the SM expectation. This helps ensure that we do not fit an artificially high value of Y t by ignoring higher-order diagrams. This represents one of the most significant uncertainties in the fit, accounting for approximately 7% of the final measurement uncertainty. It is also observed to primarily affect the lower bound of the measurement, thus reducing our ability to distinguish between values of Y t < 1.
The uncertainty in modeling the initial-and final-state radiation in the parton shower algorithm is assessed by varying the value of the renormalization scales in the initial-and final-state radiation by a factor of two. These are among the most limiting modeling uncertainties in the measurement, contributing about 8% of the total measurement uncertainty. Uncertainties for other parameters in the parton shower description are considered separately. The h damp parameter, which controls the ME to parton shower matching in POWHEG+PYTHIA, is set to the nominal value of h damp Dedicated MC samples are generated with the top quark mass varied up and down by 1 GeV from the nominal value m t = 172.5 GeV to estimate the effect of the uncertainty in the measured mass value. While this uncertainty has a significant shape effect, it ultimately accounts for only about 1% of the total measurement uncertainty. It should be noted that, although the mass and Yukawa coupling are generally treated as independent in this measurement, varying the mass will slightly modify the definition of Y t = 1. However, this effect, which is below 1%, is much smaller than the sensitivity of the measurement and can therefore be ignored.
The NNPDF sets [26] contain 100 individual variations as uncertainties. Following the approach in Ref. [13], similar variations are combined to reduce the number of variations to a more manageable set of 10 templates. The variation of the strong coupling α S used by NNPDF is treated separately from the other PDF variations. The effect of uncertainties in the PDF set is typically smaller than 1%, and together they account for 2% of the total measurement uncertainty.
The branching fraction of semileptonic b hadron decays affects the b jet response. The effect of varying this quantity within its measured precision [55] is included as an uncertainty, which has a small effect relative to other modeling uncertainties.
The momentum transfer from b quarks to b hadrons is modeled with a transfer function dependent on x b = p T (b-hadron)/p T (b-jet). To estimate the uncertainty, the transfer function is varied up and down within uncertainty of the Bowler-Lund parameter [56] in PYTHIA. The resulting effect is included by modifying event weights to reproduce the appropriate transfer function. This has a noticeable shape effect of the order 4%, but was not found to be a leading uncertainty in the fit.
In some measurements performed strictly in the context of the SM (for example, in Ref. [57]), an additional uncertainty is included to account for an observed difference in the top quark p T distribution between data and POWHEG+PYTHIA simulation. As the measurement presented here is sensitive to anomalously high values of Y t , we do not want to include any additional uncertainties which explicitly enforce agreement between SM simulation and the data, as this could reduce our sensitivity to deviations from the SM.
With this in mind, studies were performed comparing different simulations to assess whether top quark p T modelling disagreements necessitated the inclusion of any additional uncertainties. Fixed-order calculations were studied for tt production at NNLO, which generally show better agreement with the top quark p T spectrum observed in data (see, for example, Refs. [41,58]). Specifically, differential cross sections in top quark p T and M tt were studied using publicly available FASTNLO tables [59, 60], as well as multidifferential cross sections [61]. Such NNLO calculations use a different choice of dynamical scale in evaluating the top quark p T versus other kinematic variables, lending them an edge in precision over full event simulation. We find that the predictions from POWHEG+PYTHIA samples are consistent with the differential and multidifferential cross sections from Refs. [59][60][61] involving M tt and ∆y tt , within modeling uncertainties. These distributions appear consistent with the data as well, as previously observed in Ref. [41]. By comparison, the top quark p T distribution evaluated at NNLO from Refs. [59,60] shows more substantial disagreement with POWHEG+PYTHIA simulations. We conclude that the variables relevant to our measurement technique appear sufficiently well described by POWHEG+PYTHIA simulations, and differences with relevant NNLO calculations should be covered by the standard uncertainty estimation techniques. However, analyses that are more specifically sensitive to the top quark p T distribution should take care in addressing this discrepancy when using POWHEG+PYTHIA samples.

Treatment of systematic uncertainties
In this analysis, the effect of the parameter of interest Y t manifests itself as a smooth shape distortion of the kinematic distributions, as shown in Fig. 7. Although the nuisance parameters describing the sources of uncertainty should induce smooth shape effects as well, their effects are sometimes obscured by statistical noise or imprecise methods of estimation. This is noticeable for the uncertainties associated with the jet energy scale, jet energy resolution, parton shower modeling, pileup reweighting, and top quark mass. For these templates only, we apply a one-iteration LOWESS algorithm [62] to smooth the templates and remove fluctuations that may disturb the fit. The underlying-event tune and h damp uncertainties in the parton showering are small enough for their shapes to disappear into statistical noise, and are therefore treated only as normalization uncertainties.
Most templates are also symmetrized, by taking the larger effect of the up and down variations in each bin and using this magnitude for both. This step helps ensure a stable minimum in the likelihood fit, but is skipped for the templates whose natural shape effect is notably asymmetric. In the few cases where this may be an overly conservative approach, it nonetheless guarantees the performance and reliability of the minimization procedure, and has little effect on the final result.
Full or partial correlations between the 2016, 2017, and 2018 data analyses are assumed for many uncertainties. In general, the theoretically motivated uncertainties are considered fully correlated between years. Exceptions are made in cases where modeling differed between years. The PDF uncertainties cannot be correlated between 2016 and other data-taking periods, as the PDF sets used for simulation were changed to a newer version. Due to changes in the PYTHIA tune following 2016, the nominal scales used initial-state radiation and finalstate radiation differ after 2016, so those uncertainties are treated as only partially correlated between 2016 and other data-taking periods. The modeling of these uncertainties differs in the 2016 simulation, so the associated nuisance parameter in this year is either partially or fully decorrelated from those in the other years. Additionally, uncertainties whose effects disappear into statistical noise due to limited MC sample size (underlying-event tune and h damp ) are converted to uncorrelated normalization uncertainties.
Some experimental uncertainties can be broken into components, which are either fully correlated or uncorrelated between years (large jet energy scale contributions and integrated luminosity). The uncertainty in the number of pileup events is considered fully correlated as it is evaluated by varying the total inelastic cross section. For minor uncertainties from jet and lepton scale factors, which have both correlated and statistical components, a 50% correlation is assumed between years. Lastly, the jet energy resolution uncertainties are treated as uncorrelated between years.

Results
We obtain a fit result of Y t = 1.16 +0.07 −0.08 (stat) +0.17 −0.27 (syst) and an approximate upper limit at 95% CL of Y t < 1.54, where the latter is determined from the point at which −2 ln(L(Y t )) increases by an amount of 1.64 2 relative to the minimum value. For comparison, the standard model expectation based on simulated Asimov data [63] is Y t = 1 +0. 30 −0.57 (tot) with Y t < 1.47 at 95% CL. The scan of the profile likelihood test statistic used to build these intervals is shown in Fig. 8, along with a comparison to the expected behavior based on simulated Asimov data sets. We also show the agreement of data and simulation after performing the fit in Fig. 9. The minimum of the negative log likelihood occurs at a configuration with good agreement between data and simulation. The result is seen to be clearly limited by systematic uncertainties rather than statistical uncertainty. The templates for the four uncertainties with the greatest effect on the fit are shown in Fig. 10.
This result is in agreement with the previously obtained measurement in the lepton+jets final state in Ref.
[13], while obtaining a slight increase in sensitivity. Using a different decay channel and a larger data set provides a measurement complementary to the previous result.   Figure 9: The comparison between data and MC simulation at the best fit value of Y t = 1.16 after performing the likelihood maximization, with shaded bands displaying the post-fit uncertainty. The solid lines separate the three data-taking periods, while the dashed lines indicate the boundaries of the two |∆y b b | bins in each data-taking period, with M bb bin ranges displayed on the x axis. The lower panel shows the ratio of data to the simulated events in each bin, with total post-fit uncertainty bands drawn around the nominal expected bin content.

Summary
A measurement of the Higgs Yukawa coupling to the top quark is presented, based on data from proton-proton collisions collected by the CMS experiment. Data at a center-of-mass energy of 13 TeV is analyzed from the LHC Run 2, collected in 2016-18 and corresponding to an integrated luminosity of 137 fb −1 . The resulting best fit value of the top quark Yukawa coupling relative to the standard model is given by Y t = 1.16 +0.24 −0.35 . This measurement uses the effects of virtual Higgs boson exchange on tt kinematic properties to extract information about the coupling from kinematic distributions. Although the sensitivity is lower compared to constraints obtained from studying processes involving Higgs boson production in Refs.
[9] and [11], this measurement avoids dependence on other Yukawa coupling values through additional branching assumptions, making it a compelling independent measurement. This measurement also achieves a slightly higher precision than the only other Y t measurement that does not make additional branching fraction assumptions, performed in the search for production of four top quarks. The four top quark search places Y t < 1.7 at a 95% confidence level [12] while this measurement achieves an approximate result of Y t < 1.54.  Along with the intrinsic uncertainty in the EW corrections, these are the limiting uncertainties in the fit. The shaded bars represent the raw template information, while the lines show the shapes after smoothing and symmetrization procedures have been applied. In the fit, the jet energy corrections are split into 26 different components, but for brevity only the total uncertainty is shown here. Variation between years is minimal for each of these uncertainties, although they are treated separately in the fit.

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:   [4] CMS Collaboration, "Measurement of the top quark mass with lepton+jets final states using p p collisions at √ s = 13 TeV", Eur. Phys. J. C 78 (2018) 891, doi:10.1140/epjc/s10052-018-6332-9, arXiv:1805.01428.