Constraints on anomalous Higgs boson couplings to vector bosons and fermions in its production and decay using the four-lepton final state

Studies of CP violation and anomalous couplings of the Higgs boson to vector bosons and fermions are presented. The data were acquired by the CMS experiment at the LHC and correspond to an integrated luminosity of 137 fb−1 at a proton-proton collision energy of 13 TeV. The kinematic effects in the Higgs boson’s four-lepton decay H → 4l and its production in association with two jets, a vector boson, or top quarks are analyzed, using a full detector simulation and matrix element techniques to identify the production mechanisms and to increase sensitivity to the tensor structure of the Higgs boson interactions. A simultaneous measurement is performed of up to five Higgs boson couplings to electroweak vector bosons (HVV), two couplings to gluons (Hgg), and two couplings to top quarks (Htt). The CP measurement in the Htt interaction is combined with the recent measurement in the H → γγ channel. The results are presented in the framework of anomalous couplings and are also interpreted in the framework of effective field theory, including the first study of CP properties of theHtt and effectiveHgg couplings from a simultaneous analysis of the gluon fusion and top-associated processes. The results are consistent with the standard model of particle physics.


I. INTRODUCTION
The discovery by the ATLAS and CMS Collaborations of a Higgs boson (H) with a mass of ≈125 GeV [1][2][3] has confirmed the predictions of the standard model (SM) of particle physics [4][5][6][7][8][9][10]. The CMS [11][12][13][14][15][16][17][18] and ATLAS [19][20][21][22][23][24][25] experiments have set constraints on the spin-parity properties of the H boson and anomalous HVV couplings, where V stands for W, Z, and γ electroweak (EW) gauge bosons, finding its quantum numbers to be consistent with J PC ¼ 0 þþ , but leaving room for small anomalous HVV couplings. In theories beyond the SM (BSM), H boson interactions may generate several of them, which lead to new interaction tensor structures, both CP even and CP odd. The new anomalous tensor structures of the H boson interactions may also appear through loop corrections in SM processes, but the size of their contributions is beyond the current experimental sensitivity. The CP-odd anomalous couplings between the H boson and either the top quark or new BSM particles, fermions or bosons, contributing to the gluon fusion loop may generate CP violation in Hgg interactions (where g stands for gluon). Possible CP violation effects in couplings to fermions, Hff, had not been experimentally probed until recently, when the first constraints were reported by CMS [26] and ATLAS [27] in ttH production (where t stands for top quark) using the H → γγ channel.
In this paper, we study the tensor structure of the HVV, Hgg, and Htt interactions, and we search for several anomalous effects, including CP violation, using the four-lepton final state H → VV → 4l, where l ¼ μ or e. The H boson production processes considered in this paper include gluon fusion (ggH), vector boson fusion (VBF), associated production with a weak vector boson (VH, either ZH or WH), with a top quark pair (ttH), with a single top quark (tH), and with a bottom quark pair (bbH). The Feynman diagrams for these processes are shown in Figs. 1-5 and discussed in detail in Sec. II. Kinematic effects in the H boson's decay and its production in association with two jets generated in either the VBF or ggH processes, with a vector boson, or with top quarks are analyzed. Production and decay processes of the H boson are sensitive to certain anomalous contributions, or equivalently higher-dimensional operators in the effective field theory (EFT) [28], which modify the kinematic distributions of the H boson's decay products and of the particles produced in association with the H boson. Prior measurements of EW processes limit the allowed values of certain EFT operators, and the preferred EFT basis used in this paper is chosen to minimize the number of independent operators and their correlations. This allows us to reduce the number of operators to be measured. The results are also translated and reported in different, frequently used EFT bases. The effects of EFT modifications in backgrounds are found to be negligible because of the high purity of the signal peak in the four-lepton invariant mass distribution and further constraints on the backgrounds from sidebands.
Each production process of the H boson is identified using its kinematic features, and events are assigned to corresponding categories. Two categorization schemes are employed in this analysis, one designed to study Htt and Hgg and the other designed to study HVV anomalous couplings. Within each category, the matrix element likelihood approach (MELA) [29][30][31][32][33] is employed to construct observables that are optimal for the measurement of the chosen anomalous couplings, or EFT operators, including CP-sensitive observables to search for CP-violating operators. In our approach, fully simulated Monte Carlo (MC) signal samples that include anomalous couplings allow the incorporation of detector effects into the likelihood analysis, and observables explore all kinematic features of the events, including those sensitive to CP violation and to simultaneous anomalous effects in the production and decay of the H boson. These features distinguish our analysis from the recent measurements by the ATLAS Collaboration in the same H → 4l decay channel [34,35], in which case the measurements are derived from the differential distributions based on either the simplified template cross sections [28] or unfolded fiducial measurements.
We follow the formalism used in the study of anomalous couplings in the earlier analyses of CMS data in Refs. [11][12][13][14][15][16][17][18]26]. We focus on the measurements where the H boson is produced on shell; the extension to the off shell region is considered in Ref. [17], where joint constraints on the H boson width Γ H and its anomalous couplings are obtained using a partial dataset. Some of the theoretical foundations relevant to the present analysis can be found in Refs. [28][29][30][31][32][33]. The results are presented in a model-independent way, which allows interpretation in the scattering amplitude or effective Lagrangian approach, for example in the frameworks of the standard model effective field theory (SMEFT) [67][68][69][70] or pseudo observables [28].
Compared to our previous results on anomalous HVV couplings [16,17], which used a subset of the data presented here, several substantial improvements have been introduced. First, as a result of the increased number of H → 4l events, the expected 95% confidence level (CL) constraints on HVV couplings are now dominated by the tight limits from the analysis of kinematic distributions of particles produced in association with the H boson in the VBF and VH processes. Second, an improved fit implementation allows for the simultaneous measurement of up to five independent HVV, two Hgg, and two Htt couplings using the single decay channel H → 4l. The couplings are parametrized with the signal strength and the fractional cross section contributions of the anomalous couplings. A direct constraint on CP violation in the Hgg coupling is obtained for the first time by employing CP-sensitive observables. The CP violation measurement in the Htt coupling closely follows our recent measurement in the H → γγ channel [26], and the results are combined. We also perform the first study of CP properties in the Htt and effective Hgg couplings from a simultaneous analysis of the gluon fusion and top-associated processes. For all anomalous couplings, we interpret our results via the SMEFT framework in terms of HVV, Hgg, and Htt operators.
The rest of this paper is organized as follows. The phenomenology of anomalous HVV, Hgg, and Hff couplings and considerations in the EFT framework are discussed in Sec. II. The data and MC simulation, event reconstruction, and selection are discussed in Sec. III. The kinematic variables associated with the H boson's production and decay and its MELA analysis are introduced in Sec. IV. The implementation of the maximum likelihood fit is shown in Sec. V. The results are presented and discussed in Sec. VI. We conclude in Sec. VII.

II. PARAMETRIZATION OF H BOSON PRODUCTION AND DECAY PROCESSES
The goal of this study is to search for CP violation and, more generally, anomalous couplings of the H boson, in its interactions with fermions and vector bosons in the production and decay processes. These potential sources of CP violation and anomalous tensor structures of interactions may arise from BSM effects, including those considered in the EFT formulation. The dominant processes sensitive to such interactions are shown in Figs. 1-5 [28].
The main decay process considered in this paper is H → VV → 4l, with the HVV vertex shown in Fig. 1, right. The dominant H boson production mechanism is gluon fusion ggH, shown in Fig. 1, left. The dominant contribution to the gluon fusion loop comes from the top quark, with smaller contributions from the bottom quark and lighter quarks. However, contribution of new BSM states to the loop and variation of CP properties of the H boson couplings to SM quarks are also possible and are considered in this paper. The HVV vertex also appears in the vector boson fusion VBF and associated production with a weak vector boson ZH or WH, shown in Fig. 2, which are the next dominant production mechanisms of the H boson.
The production of an H boson in association with a top quark pair ttH is shown in Fig. 3, left. This is the main channel that allows to study the CP property in the H boson coupling to fermions. We also combine our results with the recent ttH measurements in the H → γγ channel [26]. The production of an H boson in association with a single top quark tH is shown in Fig. 4. This production receives contributions from both HVV and Htt couplings, but its expected cross section is smaller than that of ttH production. Both HVV and Htt couplings also contribute to the gg → ZH production mode, shown in Fig. 5. However, this gluon fusion production mode of ZH is expected to contribute only about 5% of the VH process cross section shown in Fig. 2, the dependence on anomalous HVV couplings is suppressed in this process [33], and for these reasons this production mode is neglected in this analysis. Finally, we also consider the bbH production mode shown in Fig. 3, right. However, this process does not provide kinematic features that could distinguish the CP structure    of interactions [32] or the experimental signatures that would allow its isolation from the other more dominant production mechanisms.

A. Parametrization of production and decay amplitudes
Anomalous effects in the H boson couplings to fermions, such as in the ttH and bbH production and partially in the tH and gg → ZH production, can be parametrized with the amplitude defined for each fermion type f, whereψ f and ψ f are the fermions' Dirac spinors, κ f andκ f are the corresponding coupling strengths, m f is the fermion mass, and v is the SM Higgs field vacuum expectation value. In the SM, the coupling strengths are κ f ¼ 1 andκ f ¼ 0. The presence of both CP-even κ f and CP-oddκ f couplings will lead to CP violation. In an experimental analysis of the bbH process it is not possible to resolve the κ b andκ b couplings [32], but it is possible to resolve the κ t andκ t couplings in the ttH and tH processes, which we explore in this paper. Anomalous effects in EW H boson production (VBF, ZH, and WH), ggH production, H → VV decay, and partially in the tH and gg → ZH production, are described by the HV 1 V 2 couplings. The scattering amplitude describing the interaction between a spin-zero H boson and two spin-one gauge bosons V 1 V 2 , such as ZZ, Zγ, γγ, WW, or gg, is written as where f ðiÞμν ¼ ϵ μ Vi q ν Vi − ϵ ν Vi q μ Vi ,f ðiÞ μν ¼ 1 2 ϵ μνρσ f ðiÞ;ρσ , and ϵ Vi , q Vi , and m Vi are the polarization vector, fourmomentum, and pole mass of a gauge boson i ¼ 1 or 2. The constants Λ 1 and Λ Q are the scales of BSM physics necessary to keep the κ VV i couplings unitless, and a VV 1 , a VV 2 , a VV 3 , κ VV 1 , κ VV 2 , and κ VV 3 are real numbers that modify the corresponding amplitude terms. Equation (2) describes couplings to both EW bosons and gluons, so HV 1 V 2 can stand for HVV or Hgg.
In Eq. (2), the only nonzero tree-level contributions in the SM are a ZZ 1 ≠ 0 and a WW 1 ≠ 0. In the SM, The rest of the ZZ and WW couplings are considered to be anomalous contributions, which are either small contributions arising in the SM because of loop effects or new BSM contributions. Among the anomalous contributions, considerations of symmetry and gauge invariance require κ ZZ . Therefore, there are a total of 13 independent parameters describing couplings of the H boson to EW gauge bosons and two parameters describing couplings to gluons. The presence of any of the CP-odd couplings a VV neglected. We apply additional symmetry considerations discussed below, which reduce the number of independent parameters to be measured. The chosen basis coincides with the Higgs basis [28] under SUð2Þ × Uð1Þ symmetry. The choice of the EFT operator basis is motivated by the natural parametrization in terms of observable states and, as a consequence, allows a more transparent construction of the data analysis and presentation of the results. However, we also present results in the Warsaw basis [28], using tools in Refs. [33,71] to perform the translation. Our approach with SUð2Þ × Uð1Þ symmetry is equivalent to the SMEFT formulation [28].
The parametrization in Eq. (2) is the most general one, and we apply SUð2Þ × Uð1Þ symmetry in the relationships of anomalous couplings as follows [28,33]: a WW Since m W is measured to high precision [72], we set Δm W ¼ 0, leading to a ZZ 1 ¼ a WW 1 . The latter relationship also appears under custodial symmetry [73]. Therefore, the set of 13 þ 2 independent parameters describing the HVV þ Hgg couplings can be reduced to 8 þ 2 with the above symmetry relationships.
In our measurements, we further reduce the number of independent parameters in the following way. We assume that the four loop-induced couplings a γγ 2;3 and a Zγ 2;3 are constrained to yield the SM rates of the direct decays H → γγ and Zγ. Therefore, in our analysis of EW production and H → 4l decay, we set these four couplings to zero because their allowed values are expected to have negligible effect in our coupling measurements. These four anomalous couplings have been tested in our earlier analysis of the H → VV → 4l process with run 1 data [13] and the obtained constraints were significantly looser than those from the direct decays with on shell photons. We adopt two approaches to set the relationship between the HZZ and HWW couplings. The HWW couplings do not contribute to the H → 4l decay, but they do contribute to the EW production. In this paper, the relationship between the HZZ and HWW couplings is mostly relevant for VBF production with ZZ and WW fusion. There are no kinematic differences between these two processes and it is impossible to disentangle the HZZ and HWW couplings from these data. We used approach 1 in our previous publications [16][17][18], where we set a WW i ¼ a ZZ i . Approach 2 corresponds to the SMEFT formulation with SUð2Þ × Uð1Þ symmetry.
In approach 1, we set the ZZ and WW couplings to be equal, a WW i ¼ a ZZ i . Formally, this could be considered as the relationship in Eqs. (3)- (6) in the limiting case c w ¼ 1.
As a result, we are left with four anomalous couplings to be measured, in addition to the SM-like couplings a 1 : a 2 , a 3 , κ 1 =ðΛ 1 Þ 2 , and κ Zγ 2 =ðΛ Zγ 1 Þ 2 , where we drop the ZZ superscript from the couplings. We adopt this approach both for its simplicity and to be able to relate the four anomalous couplings constrained with the H → 4l decay to the equivalent four couplings in the pseudo observables approach. This requires an independent measurement of the κ Zγ 2 term, which would otherwise be eliminated by the relationship in Eq. (7). Therefore, this approach is slightly less restrictive than the SMEFT formulation adopted in approach 2 discussed below.
In approach 2, we adopt the full set of SUð2Þ × Uð1Þ symmetry relationships in Eqs. (3)-(7) with s 2 w ¼ 0.23119 [72]. The number of independent HVV parameters is further reduced from five to four: a 1 , a 2 , a 3 , and κ 1 =ðΛ 1 Þ 2 . There is a linear one-to-one relationship between the amplitude couplings in Eq. (2) and the EFT couplings in the Higgs basis in the notation of Refs. [28,33]: Ignoring small loop-induced corrections, the above four parameters are expected to be zero in the SM. Since we set a γγ 2;3 and a Zγ 2;3 to zero, the four corresponding parameters in the EFT Higgs basis c γγ , c zγ ,c γγ , andc zγ are also zero.
In the case of Hgg couplings, the two EFT parameters are defined following the notation of Refs. [28,33] as where in the SM, c gg ¼ 0 andc gg ¼ 0, and the SM process is generated by the quark loop not accounted for in the a gg 2 coupling. Finally, in the case of the Hff couplings, the κ f andκ f parameters in Eq. (1) can be treated as EFT parameters [28], where in the SM, κ f ¼ 1 andκ f ¼ 0.

C. Parametrization of cross sections
We use the narrow-width approximation and parametrize differential cross section of the on shell H boson production process j and decay to a final state f following Refs. [28,33] as where a i are the real couplings describing the Hff, Hgg, or HVV interactions and include generically the κ i in Eqs. (1) and (2). The coefficients α ðkÞ il are in general functions of kinematic observables for the differential cross section distributions and are modeled with simulation. The total width Γ H depends on the couplings a i and potentially on the partial decay width to unobserved or invisible final states, a dependence that must be taken into account when interpreting cross section measurements in terms of couplings.
When we perform the amplitude analysis of the data, in the likelihood based on Eq. (14), we keep all cross terms in the expansion of powers of Λ −N with N ¼ 0, 2, 4, 6, 8, where formally each dimension-6 operator, or anomalous a i coupling, receives the Λ −2 contribution in Eq. (2), even if not explicitly shown, while the SM tree-level coupling carries no such contribution. While this may create inconsistency in the EFT approach with the Λ −N terms kept or neglected from the higher-dimension contributions, this allows us to keep the likelihood positive definite, which is an important consideration in the experimental analysis of the data discussed in Sec. V. Because interference contributions may become negative, dropping certain terms in the expansion may lead to negative probability. The importance of N ¼ 4, 6, 8 contributions may be considered as testing whether the current precision is sufficient to treat our results within the EFT approach, and we leave this test to the interpretation of the results. However, regardless of EFT validity, our results are presented in a fully selfconsistent formulation of amplitude decomposition, which can either be translated to the EFT interpretation or treated as a test of consistency of the data with the SM, including the search for new sources of CP violation.

D. Parametrization of the signal strength and cross section fractions
We present the primary results in terms of cross sections, or equivalently, signal strengths μ j ¼ σ j =σ SM j , and the fractional contributions f ai of the couplings a i to cross sections ð P mn α mn a m a n Þ of a given decay process. The ratios of couplings entering Eq. (14) can be expressed through f ai , and the common factors, such as the total width Γ H and the SM-like coupling squared, are absorbed into the signal strength. This formulation with μ j and f ai allows the presentation of experimental results in the most direct way, with a minimal and complete set of parameters describing the given processes. This approach has several convenient features. The cross sections and their ratios are invariant with respect to the coupling convention, such as the scaling in Eqs. (8)- (13). The cross section fractions f ai reflect kinematic features in either production or decay in a direct way. They are conveniently bound between −1 and þ1, and most systematic uncertainties cancel in the ratio.
The cross section fraction for Hff couplings is defined as Similarly, the cross section fraction for Hgg couplings is defined as Both definitions incorporate the relative sign of the possible BSM CP-odd and SM-like CP-even couplings. They are based on the observation that the cross sections of the H → gg process are equal for a gg 2 ¼ 1 and a gg 3 ¼ 1, as are the cross sections of the H → ff process for κ f ¼ 1 and κ f ¼ 1 in the limit of m f ≪ m H . We note that f ggH a3 is defined following the convention that a gg 2 and a gg 3 absorb both pointlike interactions and quark contributions to the loop. Following Ref.
[33], the f ggH a3 measurement can also be interpreted in terms of f Hff CP under the assumption that only the top and bottom quarks contribute to gluon fusion where the signs of f Hff CP and f ggH a3 are equal, and α Hff is an effective parameter sometimes used to describe the CP-odd contribution to the H boson Yukawa couplings. A more detailed analysis of the gluon fusion loop could be performed without the assumption that only the top and bottom quarks contribute.
The cross section fractions in the HVV couplings of the H boson to EW gauge bosons require more parameters. Since in both of our approaches the HWW couplings are expressed through other a VV i couplings following Eqs. (3)- (6), and because we prefer that our definitions not depend on parton distribution functions (PDFs) and other effects that involve measurement uncertainties, we use the H → ZZ=Zγ Ã =γ Ã γ Ã → 2e2μ decay process to define the cross section fractions as where the α ð2e2μÞ ii coefficients are introduced in Eq. (14). The numerical values of these coefficients are given in Table I, where they are normalized with respect to the  α   ð2e2μÞ  11 coefficient, corresponding to the cross section calculated for a 1 ¼ 1. The α ð2e2μÞ ii are the cross sections for a VV i ¼ 1, which are different in the two approaches of the coupling relationship as a result of Eq. (7) adopted in approach 2. The cross section fractions in Eq. (18) can be converted to coupling ratios as The measured values of μ j and f ai should be sufficient to adopt them in the fits for EFT parameters jointly with the data from other H boson, top quark, and EW measurements. They allow constraints on the κ i and a i couplings in Eqs. (1) and (2). However, it is required to perform a simultaneous measurement of all production and decay channels of the H boson, including unobserved and invisible channels, as they contribute to the total width in Eq. (14). In this paper, we present only a limited interpretation of our data in terms of couplings by making certain assumptions about their relationship. We leave more extensive interpretation to a future combination with other channels.

III. THE CMS DETECTOR, DATASETS, AND EVENT RECONSTRUCTION
The H → 4l decay candidates are produced in protonproton (pp) collisions at the LHC and are collected and reconstructed in the CMS detector [74]. The data sample used in this analysis corresponds to integrated luminosities of 35.9 fb −1 collected in 2016, 41.5 fb −1 collected in 2017, and 59.7 fb −1 collected in 2018, for a total of 137 fb −1 collected during Run 2 at a pp center-of-mass energy of 13 TeV.
The CMS detector comprises a silicon pixel and strip tracker, a lead tungstate crystal electromagnetic calorimeter (ECAL), and a brass/scintillator hadron calorimeter, each composed of a barrel and two end cap sections, all within a superconducting solenoid of 6 m internal diameter, providing a magnetic field of 3.8 T. Extensive forward calorimetry complements the coverage provided by the barrel and end cap detectors. Outside the solenoid are the gas-ionization detectors for muon measurements, which are embedded in the steel flux-return yoke. A detailed description of the CMS detector can be found in Ref. [74].
Events of interest are selected using a two-tiered trigger system. 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 fixed latency of about 4 μs [75]. The second level, known as the high-level trigger, 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 [76].

A. Event reconstruction and selection
The selection of 4l events and associated particles closely follows the methods used in the analyses of the run 1 [12,13] and run 2 [16,17,77,78] datasets. The main triggers for the run 2 analysis select either a pair of electrons or muons, or an electron and a muon, passing loose identification and isolation requirements. The transverse momentum (p T ) for the leading electron or muon is required to be larger than 23 or 17 GeV, while that of the subleading lepton is required to be larger than 12 or 8 GeV, respectively. To maximize the signal acceptance, triggers requiring three leptons with lower p T thresholds and no isolation requirement are also used, as well as isolated single-electron and single-muon triggers with thresholds of 27 and 22 GeV in 2016, or 35 and 27 GeV in 2017 and  (18), and the translation coefficients α ii =α 11 in this definition with the relationship a ZZ i ¼ a WW i (approach 1), and with the SMEFT relationship according to Eqs. (3)-(7) (approach 2). In the case of the κ 1 and κ Zγ 2 couplings, the numerical values Λ 1 ¼ Λ Zγ 1 ¼ 100 GeV are adopted in this calculation to make the coefficients have the same order of magnitude and the negative sign indicates the convention in Eq. (18) adopted earlier [13]. In approach 2, κ Zγ 2 is a dependent parameter expressed through Eq. (7) and does not require a translation coefficient.
---2018, respectively. The overall trigger efficiency for simulated signal events that pass the full selection chain of this analysis is larger than 99%. Event reconstruction is based on the particle-flow algorithm [79], which exploits information from all the CMS subdetectors to identify and reconstruct individual particles in the event. The particle-flow candidates are classified as charged or neutral hadrons, photons, electrons, or muons, and they are then used to build higher-level objects, such as jets, and to calculate the lepton isolation quantities. Electrons or muons are reconstructed within the geometrical acceptance defined by a requirement on the pseudorapidity jηj < 2.5 or 2.4 and p T > 7 or 5 GeV, with an algorithm that combines information from the tracker and the ECAL or muon system, respectively. Muons are selected from a list of reconstructed muon track candidates by applying minimal requirements on the track in both the muon system and inner tracker system, and taking into account compatibility with small energy deposits in the calorimeters.
To discriminate between leptons from prompt particle decays and those arising from hadron decays within jets, an isolation variable is calculated for electrons and muons [78]. An isolation requirement is imposed on the muons. Electrons are identified using a multivariate discriminant which includes observables sensitive to the presence of bremsstrahlung along the electron trajectory, the geometrical and momentum-energy matching between the electron trajectory and the associated cluster in the ECAL, the shape of the electromagnetic shower in the ECAL, and variables that discriminate against electrons originating from photon conversions. This discriminant also includes the isolation to suppress electrons originating from EW decays of hadrons within jets [78]. A dedicated algorithm is used to collect the final-state radiation of leptons [77].
The jets are clustered using the anti-k T jet finding algorithm [80,81] with a distance parameter of 0.4. The jet momentum is determined as the vector sum of all particle momenta in the jet. Jets must satisfy p T > 30 GeV and jηj < 4.7 and must be separated from all selected lepton candidates and any selected final-state radiation photons with a requirement on the parameter ΔRðl=γ; jetÞ > 0.4, where ðΔRÞ 2 ¼ ðΔηÞ 2 þ ðΔϕÞ 2 . Jets are b tagged using the DeepCSV algorithm [82], which combines information about impact parameter significance, the secondary vertex, and jet kinematics.
The reconstructed vertex with the largest value of summed physics-object p 2 T is taken to be the primary pp interaction vertex. The physics objects are the jets, clustered using the jet finding algorithm [80,81] with the tracks assigned to candidate vertices as inputs, and the associated missing transverse momentum, taken as the negative vector p T sum of those jets. In order to suppress muons from in-flight decays of hadrons and electrons from photon conversions, leptons are rejected if the ratio of their impact parameter in three dimensions, computed with respect to the primary vertex position, to their uncertainty is greater than four.
We consider three mutually exclusive lepton flavor channels: H → VV → 4e, 4μ, and 2e2μ. At least two leptons are required to have p T > 10 GeV, and at least one is required to have p T > 20 GeV. All four pairs of oppositely charged leptons that can be built with the four leptons are required to satisfy m l þ l − > 4 GeV, regardless of lepton flavor, to further suppress events with leptons originating from hadron decays in jet fragmentation or from the decay of low-mass resonances. The V candidates are formed with pairs of leptons of the same flavor and opposite charge that pass the requirement 12 < m l þ l − < 120 GeV, where m 1 is the invariant mass of the V candidate that is closest to the nominal Z boson mass and m 2 is the mass of the other one. A value of m 1 > 40 GeV is required. The reconstructed four-lepton invariant mass, m 4l , distribution in the region between 70 and 170 GeV is shown in Fig. 6. The m 4l region between 105 and 140 GeV is considered in this analysis, which is wide enough to use sidebands for constraining the background normalization in the later fitting procedure [78].

B. Event categorization
In order to perform a dedicated study of a particular kinematic topology, events are further split into several mutually exclusive categories based on the presence of other particles produced in association with the H boson candidate [78]. Two independent categorization, schemes 1 and 2, discussed below, are followed in this study. Scheme 1 targets Htt and Hgg anomalous couplings, while scheme 2 targets HVV anomalous couplings. We use the values of kinematic discriminants and other selection requirements to perform the categorization. The definition of these discriminants can be found in Refs. [16,17,77,78] and is further discussed in Sec. IV. They are calculated using the MELA while employing the matrix elements at leading order (LO) in quantum chromodynamics (QCD). These discriminants use full kinematic information from the H boson and from associated jet production and are labeled to indicate a specific topology (1jet, 2jet) and production mechanism (VBF, WH, ZH), which is discriminated against the dominant gluon fusion process: D VBF 1 jet , D VBF 2 jet , D ZH 2jet , and D WH 2jet . The D 2jet discriminants are calculated using both SM and anomalous coupling hypotheses, leading to a set D i 2jet , all of which are tested in order to maintain high efficiency of VBF and VH categorization in the presence of anomalous couplings. The discriminants defined for the two-jet topology are illustrated in Fig. 7, where the expected distributions are based on the MC signal simulation discussed in Sec. III C and the background estimate in Sec. III D. To enhance the signal to background ratio in this illustration in Fig. 7, a selection of D bkg > 0.7 is applied. This observable uses information from the lepton kinematic distributions and does not use information from associated jets, as also discussed in Sec. IV.
In categorization scheme 1, the Htt and Hgg anomalous couplings are targeted. The categories and selection criteria are identical to the first step of the categorization scheme in Ref. [78] and are optimized to measure the rates of H boson production modes. Because anomalous Htt and Hgg couplings have only a small effect on the fractions of ttH and ggH events in each category, the optimization based on SM kinematic distributions used for the study in Ref. [78] remains optimal here. The sequential selection criteria in scheme 1 are as follows: (i) The VBF-2jet category requires exactly four leptons. In addition, there must be either two or three jets of which at most one is b tagged, or at least four jets and no b-tagged jets. Finally, D VBF 2jet > 0.5 is required.
(ii) The VH-hadronic category requires exactly four leptons. In addition, there must be either two or three jets, or at least four jets and no b-tagged jets.
The VH-leptonic category requires no more than three jets and no b-tagged jets in the event, and exactly one additional lepton or one additional pair of opposite-sign same-flavor leptons. This category also includes events with no jets and at least one additional lepton. (iv) The ttH-hadronic category requires at least four jets, if one is a b-tagged jet, and no additional leptons. (v) The ttH-leptonic category requires at least one additional lepton in the event. (vi) The VBF-1jet category requires exactly four leptons, exactly one jet and D VBF 1jet > 0.7. (vii) The untagged category consists of the remaining events. The number of events expected from signal simulation and background estimation are shown along with the observed number of events for each scheme 1 category in Table II. The distributions of observed events (data points) and expectation (histograms) for maxðD VBF;i 2jet Þ (middle) and max ðD WH;i 2jet ; D ZH;i 2jet Þ (right), where maximum is evaluated over the SM and the four anomalous coupling hypotheses i, described in the legend (left). Only events with at least two reconstructed jets are shown, and the requirement D bkg > 0.7 is applied in order to enhance the signal contribution over the background, where D bkg is calculated using decay information only. The expectation is shown for the total distribution, including background and all production mechanisms of the H boson, and for the VBF (middle) and VH (right) signals, which are enhanced in the region above 0.5, indicated with the vertical dashed line.
In categorization scheme 2, which targets anomalous HVV couplings, the categorization sequence is modified in three ways in order to be more sensitive to the HVV couplings. First, the D 2jet discriminants calculated using the SM hypothesis for VBF or VH are less sensitive to anomalous VBF or VH production, so the selection for the VBF-2jet and VH-hadronic categories is modified to be more efficient for BSM hypotheses. Second, the ttH categories are dropped and these events are merged into the untagged category; ttH forms a small background to VBF and VH. Third, a boosted category is added. This category, adapted from the second (and finer) categorization scheme in Ref. [78], is designed for events where not all associated particles are fully reconstructed, so that the full kinematic information cannot be used to measure anomalous couplings. After these modifications, scheme 2 contains six categories, with sequential selection criteria as follows: (i) The VBF-2jet category requires exactly four leptons.
In addition, there must be either two or three jets of which at most one is b tagged, or at least four jets and no b-tagged jets. Finally, maxðD VBF;i 2jet Þ > 0.5 using either the SM or any of the four BSM signal hypotheses (i) for the VBF production is required. See Fig. 7 (middle) for illustration. (ii) The VH-hadronic category requires exactly four leptons. In addition, there must be either two or three jets, or at least four jets and no b-tagged jets. Finally, max ðD WH;i 2jet ; D ZH;i 2jet Þ > 0.5 using either the SM or any of the four BSM signal hypotheses (i) for the VH production is required. See Fig. 7 (right) for illustration. (iii) The VH-leptonic category requires no more than three jets and no b-tagged jets in the event, and exactly one additional lepton or one additional pair of opposite-sign same-flavor leptons. This category also includes events with no jets and at least one additional lepton. (iv) The VBF-1jet category requires exactly four leptons, exactly one jet and D VBF 1jet > 0.7. (v) The boosted category requires exactly four leptons, three or fewer jets, or at least four jets and no b-tagged jets, and the transverse momentum of the four-lepton system p 4l T > 120 GeV. (vi) The untagged category consists of the remaining events. The number of events expected from signal simulation and background estimation are shown along with the observed number of events for each scheme 2 category in Table III. The events in each category in either schemes 1 or 2 are further characterized with several observables using the kinematic features of the H boson decay and associated particles, as discussed in Sec. IV.

C. Monte Carlo simulation
Monte Carlo simulation is used to model signal processes, which involve the H boson, and background processes in pp interactions at the LHC and their reconstruction in the CMS detector. All MC samples are interfaced with PYTHIA8 [83] to simulate parton showering and multiparton interactions, using version 8.230 for all years with the CUETP8M1 tune [84] for the simulation of the 2016 datataking period, and the CP5 tune [85] for the simulation of the 2017 and 2018 data taking periods. The NNPDF 3.0 parton distribution functions are used [86]. Simulated events include the contribution from additional pp interactions within the same or adjacent bunch crossings (pileup) and are weighted to reproduce the observed pileup TABLE II. The numbers of events expected in the SM for different H signal (sig) and background (bkg) contributions and the observed number of events in each category defined in scheme 1 targeting Hff and Hgg anomalous couplings. The ttH signal expectation is quoted for the SM and anomalous coupling (κ t ¼ 0,κ t ¼ 1.6) scenario, both generated with the same cross section. distribution in data. The MC samples are further processed through a dedicated simulation of the CMS detector based on GEANT4 [87]. The JHUGen7.3.0 [29-33] MC program is used to simulate all anomalous couplings in the H boson production and H → ZZ=Zγ Ã =γ Ã γ Ã → 4l decay as discussed in Sec. II. The MELA [29-33] package contains a library of matrix elements from JHUGen for the signal and MCFM7.0.1 [88] for the background; these matrix elements are used to apply weights to events in any MC sample to model any other set of anomalous or SM couplings.
The SM production of the H boson through VBFVBF, in association with a W or Z boson, or with a t t pair is simulated using both JHUGen at LO in QCD and POWHEG2 [89][90][91][92][93] at next-to-leading order (NLO) in QCD. Production in association with a bb pair or single top quark is simulated at LO in QCD via JHUGen. In the VBF, VH, and ttH production modes, the JHUGen and POWHEG simulations are explicitly compared after parton showering in the SM case, and no significant differences are found in kinematic observables. Therefore, the JHUGen simulation is adopted to describe kinematic distributions in the VBF, VH, ttH, tH, and bbH production modes with anomalous couplings, with the expected yields scaled to match the SM theoretical predictions [28] for inclusive cross sections and POWHEG simulation for categorization of events based on associated particles in the SM. There are no observable anomalous effects in kinematic distributions of the bbH process [32], but we keep this process in modeling its event contribution. The considered VH process does not include gg → ZH production, which is expected to contribute about 5% of the VH cross section and is therefore neglected in this analysis. This process has been studied with JHUGen, including anomalous HVV and Hff couplings, and it was found that the dependence on anomalous HVV couplings is suppressed [33].
Gluon fusion production is simulated with the POWHEG2 event generator at NLO in QCD. The kinematic features of events produced in gluon fusion with two associated jets are also modified by anomalous Hgg couplings. These effects are studied using MADGRAPH5_aMC@NLO 2.6.0 [63,94] and JHUGen. Simulation with the MINLO [95] program at NLO in QCD is used for evaluation of systematic uncertainties related to the modeling of two associated jets. The relationship between the Hff and Hgg couplings follows JHUGen with the relative sign of CP-odd and CP-even coefficients opposite to that assumed in TABLE III. The numbers of events expected in the SM for different H signal (sig) and background (bkg) contributions and the observed number of events in each category defined in scheme 2 targeting HVV anomalous couplings. The EW (VBF, WH, and ZH) signal expectation is quoted for the SM and four anomalous coupling (a 3 =a 2 =κ 1 =κ Zγ 2 ) scenarios f ai ¼ 1, all generated with the same total EW production cross section.  MADGRAPH5_aMC@NLO2.6.0, as discussed in Ref. [33]. The sign convention of the photon field in JHUGen7.3.0 is opposite to that in MADGRAPH5_aMC@NLO2.6.0, which leads to the opposite sign of the HZγ couplings. This sign convention depends on the sign in front of the gauge fields in the covariant derivative. In all of the above cases, the subsequent decay H → ZZ=Zγ Ã =γ Ã γ Ã → 4l is modeled with JHUGen. All signal processes have been generated under the assumption that the H boson mass is m H ¼ 125 GeV. This value has been used in calculations in Secs. II and III. However, in the analysis of the data discussed in Secs. IV and VI, the m H ¼ 125.38 GeV value from Ref. [96] is used. The m 4l parametrization, the cross sections and the branching fractions of all processes [28] are adjusted according to m H ¼ 125. 38 GeV, but the effect on other kinematic distributions of the H boson decay products and associated particles is neglected owing to the small difference between the two m H values.

D. Background modeling
The main background in this analysis, qq → ZZ=Zγ Ã = γ Ã γ Ã → 4l, is estimated from NLO simulation with POWHEG. A fully differential cross section has been computed at next-to-next-to-leading order (NNLO) in QCD [97] and the NNLO/NLO QCD correction is applied as a function of m 4l . The gg → ZZ=Zγ Ã =γ Ã γ Ã → 4l background process is simulated with MCFM7.0.1 [88,[98][99][100] at LO in QCD. The cross section of this background process is corrected with an NNLO K factor as a function of m 4l [101][102][103], assuming that the signal and background processes have the same correction for higher orders in QCD. The EW background includes the vector boson scattering and VZZ processes, generated within the JHUGen framework by adopting the MCFM matrix elements for the background processes. The EW background also incorporates other VVV, ttVV, and ttV processes, which are generated with MADGRAPH5_aMC@NLO. Other background contributions are estimated using control samples in reconstructed data without relying on simulation. Different sources of leptons such as leptons originating from decays of heavy flavor quarks or light mesons may produce additional background to the H boson signal. We denote this background collectively as the Z þ X background and employ a data-driven method for its estimation. The same method has been used in the analyses of run 1 [12,13] and run 2 [16,17,77,78] datasets. The lepton misidentification rates are first derived using Z þ 1l control regions with relaxed selection requirements on the third lepton, and the extracted rates are then applied in Z þ 2l control regions, where the two additional leptons with relaxed selection requirements have the same lepton flavor of equal or opposite charge [78].

IV. KINEMATIC EFFECTS IN PRODUCTION AND DECAY OF THE H BOSON
Kinematic distributions of particles produced in the H boson decay or in association with the H boson production are sensitive to the quantum numbers and anomalous couplings of the H boson. Four main production topologies are studied: ggH, VBF, ZH, or WH, and ttH or tH, as illustrated in Fig. 8.
In the H → VV → 4l decay, shown in Fig. 8, lower right, eight observables fully characterize the kinematic distributions of the decay products and the orientation of the decay frame with respect to the production axis . Sets of observables Ω prod for the ggH, VBF, VH, and ttH production processes are defined in a similar way [31,32], as shown in Fig. 8. As a result, 13 or more kinematic observables can be defined for the associated production process, with subsequent H boson decay to a four-fermion final state. The MELA is designed to reduce the number of observables to a minimum, while retaining all essential information.

A. Kinematic discriminants
Full kinematic information from each event, using either the H → VV → 4l decay or associated particles in its production, is extracted using discriminants from matrix element calculations using the MELA package. The discriminants use a complete set of mass and angular input observables Ω [29,31,32] to describe kinematic distributions at LO in QCD. Full reconstruction of the four-lepton decay chain and associated particles is employed in the matrix element calculation, following the selection chain discussed in Sec. III. Events with partial reconstruction of associated particles are retained in the analysis by using other kinematic observables, such as the transverse momentum of the reconstructed H boson. In the case of the ttH topology, where full reconstruction of the full decay chain of the top quarks is a challenging task, an approximation to the matrix element approach is achieved with machine learning [33].
Two types of discriminants are defined for either the production process, the decay process, or the full production þ decay process. These discriminants are where the probability density P of a certain process is calculated using the full kinematic description characterized by Ω for the processes denoted as "sig" for a signal model and "alt" for an alternative model, which could be an alternative H boson production mechanism (used to categorize events), background (to isolate signal), or an alternative H boson coupling model (to measure coupling parameters). The "int" label refers to the interference between the two model contributions. The probability densities P are calculated from the matrix elements provided by the MELA package and are normalized to give the same integrated cross section for both processes in the relevant phase space. This normalization leads to a balanced distribution of events in the range between 0 and 1 for the D alt discriminants, or between −1 and 1 for D int .
In the special case where the D int is calculated between CP-even and CP-odd models, it is denoted as D CP . The D CP observable is CP odd, and a forward-backward asymmetry in its distribution would indicate CP violation. This motivates the index "CP." When events are split into the VBF-1/2jet and VHhadronic categories, a set of discriminants D 1=2jet is constructed, following Eq. (20), where P sig corresponds to the signal probability density for the VBF (WH or ZH) production hypothesis in the VBF-tagged (VH-tagged) category, and P alt corresponds to that of H boson production in association with two jets via gluon fusion. When more than two jets pass the selection criteria, the two jets with the highest p T are chosen for the matrix element calculations. Thereby, the D 1=2jet discriminants separate the target production mode of each category from gluon fusion production, in all cases using only the kinematic properties of the H boson and two associated jets. The application of the D 1=2jet discriminants is described in Sec. III, where we introduce four types of discriminants D VBF 1jet , D VBF;i 2jet , D ZH;i 2jet , and D WH;i 2jet , with the SM and the four anomalous coupling hypotheses i considered in the signal model.
Several arrays of observables ⃗ x are defined in each category of events, uniquely targeting kinematic features of each category, and are listed in Table IV. One observable, D bkg , is common to most production categories in both schemes 1 and 2. This observable is calculated using Eq. (20) and is designed to separate signal from the dominant background production of four leptons. The P alt probability density is calculated for the dominant qq → 4l background process. The signal and background probability densities include both the matrix element probability based on the four-lepton kinematic properties from MELA and the empirical m 4l probability density parametrization extracted from the simulation of detector effects. In the VBF-2jet and VH-hadronic categories in scheme 2, the observable D EW bkg is a modified version of D bkg which includes the jet information. In this case, P sig and P alt still include the m 4l probability parametrization and four-lepton kinematic information, but they also include kinematic information for the two associated jets. The P alt probability density represents the EW and QCD background processes 4l þ 2jets, while P sig represents the EW H production processes summed together, VBF, WH, and ZH. The D bkg or D EW bkg calculation employs the SM hypothesis for the signal, while BSM kinematic information is incorporated in the observables discussed next.

B. Observables targeting anomalous Htt
and Hgg couplings In scheme 1, designed to study anomalous Htt and Hgg couplings, seven event categories are used. In the untagged, VBF-1jet, VH-leptonic, and VH-hadronic categories, only one observable D bkg is used. These categories do not provide additional information for separating CP-even and CP-odd contributions in the Htt and Hgg couplings, but are included in the fit in order to constrain the rates of the processes. The probability density parametrization of D bkg in these categories is not sensitive to the CP structure of either Hff or HVV interactions.
There is rich kinematic information in ttH production because of the sequential decay of the top quarks, as discussed further in Ref. [32]. While it is possible to construct observables, as defined in Eqs. (20) and (21), with matrix element techniques [32], we adopt a machine learning approach to account for partial reconstruction and possible permutations of the jets. A boosted decision tree (BDT) classifier is trained to separate CP-even, corresponding to the κ t coupling, and CP-odd, corresponding to theκ t coupling, contributions independently in the ttHleptonic and ttH-hadronic categories. The discriminant D ttH 0− is obtained with this approach as the best approximation to Eq. (20), provided the full kinematic information is made available in the calculation [32]. This technique still ensures that the maximal information is retained in the discriminant and is based on the same matrix element used in simulation.
We achieve the full kinematic information in the D ttH 0− calculation by including the following observables in the training: the four-momenta of the reconstructed H boson and of the six jets with the largest p T , as well as the b-tagging scores of the six jets for resolving their permutation. In addition, in the ttH-leptonic category, the lepton multiplicity and the four-momentum of the highest-p T lepton not originating from the H boson decay are used as input to the BDT classifier. It is not possible to construct the D ttH CP discriminant, corresponding to Eq. (21), without tagging the flavors of the jets, including distinguishing IV. The list of kinematic observables used for category selection and fitting in categorization schemes 1 and 2. Only the main features involving the kinematic discriminants in the category selection are listed, while complete details are given in Sec. III. The untagged category includes the events not selected in the other categories.

Category Selection
Observables ⃗ x for fitting CP quarks from antiquarks [32]. Alternatively, in the leptonic decay of both top quarks one could use the charges of the leptons, but the efficiency of such a method is very low. Therefore, in the two ttH categories, two observables are used: The distributions of these two discriminants are shown in Fig. 9.
The analysis of gluon fusion production with associated jets is performed in the VBF-2jet category. There are two discriminants that are sensitive to CP-even terms, corresponding to the a gg 2 coupling, and to CP-odd terms, corresponding to the a gg 3 coupling, D ggH 0− and D ggH CP , following Eqs. (20) and (21), respectively. The matrix element for the gluon fusion H boson production in association with two jets includes three possible initial states: quark-quark, quark-gluon, and gluon-gluon. Only the quark-quark initial state is used to calculate these discriminants, because this configuration corresponds to the gluon scattering topology sensitive to CP properties of the Hgg coupling [32], by analogy with the weak vector boson scattering process. The jets in the other configurations, quark-gluon and gluon-gluon, are more likely to be initiated from gluon radiation or splitting and are less likely to carry information about the CP properties. For similar reasons, we include the D VBF 2jet discriminant as one of the observables. This discriminant allows us to isolate the VBF-like topology of the events, which is more characteristic of the quark-quark-initiated process.
As a result, in the VBF-2jet category in scheme 1 the observables ⃗ x ¼ fD bkg ; D VBF 2jet ; D ggH 0− ; D ggH CP g are used, as summarized in Table IV. The distributions of the three observables D bkg ; D ggH 0− , and D ggH CP in this category are shown in Fig. 10, while the D VBF 2jet observable is shown in Fig. 7. It has been shown [33,44] that the azimuthal angle between the two jets carries information similar to D ggH 0− and D ggH CP , but the latter two are better in terms of  performance and practical application to parametrization in the fit discussed in Sec. V.

C. Observables targeting anomalous HVV couplings
In scheme 2, designed to study anomalous HVV couplings, six event categories are used. Two of these categories, VBF 2jet and VH hadronic, target full reconstruction of the associated jets in EW production of the H boson. Therefore, the full matrix element calculation using both production and decay information is employed, as discussed below and summarized in Table IV. Three other categories, boosted, VBF 1jet, and VH leptonic, also target EW production, but without full reconstruction of the associated particles. Therefore, in these categories matrix element calculations are not employed, and instead the transverse momentum of the H boson candidate p 4l T is used as the second observable. Anomalous couplings in EW production lead to a harder p 4l T spectrum. Finally, in the untagged category, dominated by the ggH events without two associated jets, matrix element calculations using H → VV → 4l decay information are employed.

Λ1
, respectively, where the index "dec" indicates that only the four-lepton decay information is used. Among the latter four interference discriminants, it was found that the two discriminants corresponding to the κ 1 =ðΛ 1 Þ 2 and κ Zγ 2 =ðΛ Zγ 1 Þ 2 couplings are strongly correlated with D dec Λ1 and D Zγ;dec

Λ1
, and therefore these two interference discriminants are not used. This observation follows from the fact that the a 1 and κ 1 =ðΛ 1 Þ 2 or κ Zγ 2 =ðΛ Zγ 1 Þ 2 couplings correspond to the same tensor structure in Eq. (2) and differ only in the q 2 V dependence. The remaining two interference discriminants, D dec int and D dec CP , corresponding to the a 2 and a 3 alternative couplings, respectively, are employed in the fit.
In the VBF-2jet and VH-hadronic categories, the system of six discriminants discussed above is extended to include both production and decay information, because these categories allow full reconstruction of associated particles. The same four types of discriminants of the untagged category following Eq. (20) are used, namely D VBFþdec in the VH-hadronic category. Here the index "VBF+dec" or "VH+dec" indicates that both production and decay information is used, which means that the kinematic information from the associated jets and the four leptons are utilized in the VBF or VH matrix element calculations. In the case of the VH process, the matrix elements of the WH and ZH processes are summed. There are more interference discriminants in cases where anomalous couplings appear both in production and decay. However, using interference discriminants with production information only following Eq. (21) is the better approach if one has to limit the number of discriminants. Therefore, D VBF int and D VBF CP are used in the VBF-2jet category, and D VH int and D VH CP are used in the VH-hadronic category. Distributions of events for the observables ⃗ x in scheme 2 are illustrated in Figs. 11-13. Here and in Figs. 9 and 10 the expected distributions are based on signal MC simulation discussed in Sec. III C and the background estimate in Sec. III D, where cross sections of all processes, including those with the BSM couplings, are set to the SM expectations. The full list of kinematic observables employed in the fit in each category is summarized in the third column of Table IV.

V. IMPLEMENTATION OF FITTING AND ASSOCIATED UNCERTAINTIES
In the analysis of the Hff, Hgg, or HVV anomalous couplings, the events are split into a total of 63 (in scheme 1, designed to study Hff and Hgg) or 54 (in scheme 2, designed to study HVV) categories according to the seven or six production categories, three lepton flavor combinations (4e, 4μ, and 2e2μ), and three data periods (2016, 2017, and 2018). Each event is characterized by its discrete category k and set of input observables ⃗ x, as discussed in detail in Sec. IV and summarized in Table IV. The observed distributions of events across the discriminants ⃗ x are illustrated in Figs. 9-13 and are compared to the expected distributions in SM and BSM. However, quantitative characterization of these distributions requires a careful analysis of the multidimensional space of observables ⃗ x and categories k, which is discussed below. Preparation of this analysis was performed in a blind way, which means that observed distributions of events were not examined until all details of the fit discussed below were finalized.

A. Likelihood parametrization
We perform an extended maximum likelihood fit [104] in which the probability density is normalized to the total event yield in each category k as a sum over all signal processes j and background processes i according to where μ j is the ratio of the observed cross section to the SM expectation, ⃗ f j is the set of unconstrained parameters describing kinematic distributions in a given process, defined in Eqs. (15), (16), (18), and ⃗ ξ jk are the constrained nuisance parameters reflecting the uncertainties in the above parametrization.
In the case of the gluon fusion process, μ j ¼ μ ggH and ⃗ f j ¼ f ggH a3 . The dependence on the CP-sensitive parameter f ggH a3 appears only in the VBF-2jet category where correlation of the two associated jets is explored. Therefore, in this category a dedicated simulation of the H boson production with two associated jets with MADGRAPH5_aMC@NLO is used, and checked against the JHUGen and MINLO simulation. The interpretation of this process in terms of fermion couplings appearing in the (middle middle), D dec CP , and D dec int . The last two plots are shown in the boosted category: D bkg (lower middle) and p 4l T with the requirement D bkg > 0.7 and overflow events included in the last bin (lower right). Observed data, background expectation, and five signal models are shown in the plots, as indicated in the legend in Fig. 7 (left). In several cases, a sixth signal model with a mixture of the SM and BSM couplings is shown and is indicated in the legend explicitly. gluon fusion loop is discussed in Sec. VI C. In cases where the SM fermions are assumed to dominate the gluon fusion loop, the μ ggH and f ggH a3 parameters are correlated to μ ttH and f Htt CP in the ttH process through Eq. (17). A more general case, when both SM fermions and heavy BSM particles contribute to the loop, is also considered. In all cases the relationship between the Hff and Hgg couplings follows JHUGen with the relative sign of CPodd and CP-even coefficients opposite to that assumed in MADGRAPH5_aMC@NLO, as discussed in Ref. [33]. (middle middle), D VBF CP , and D VBF int . The last two plots are shown in the VBF-1jet category: D bkg (lower middle) and p 4l T with the requirement D bkg > 0.7 and overflow events included in the last bin (lower right). Observed data, background expectation, and five signal models are shown in the plots, as indicated in the legend in Fig. 7 (left). In several cases, a sixth signal model with a mixture of the SM and BSM couplings is shown and is indicated in the legend explicitly.
In the case of the ttH process, μ j ¼ μ ttH and ⃗ f j ¼ f Htt CP . The tH production is expected to contribute about 5% of the total ttH signal in this analysis, as shown in Table II, and, therefore, the exact treatment of this process is currently not important. This production mode depends on both HVV and Hff couplings. The anticipated small yield may be somewhat larger if the expected destructive interference between the HVV and Hff couplings in the tH process becomes constructive owing to a modification of these couplings. In this analysis, we also introduce a (middle middle); D VH CP , and D VH int . The last two plots are shown in the VH-leptonic category: D bkg (lower middle) and p 4l T with the requirement D bkg > 0.7 and overflow events included in the last bin (lower right). Observed data, background expectation, and five signal models are shown in the plots as indicated in the legend in Fig. 7 (left). In several cases, a sixth signal model with a mixture of the SM and BSM couplings is shown and is indicated in the legend explicitly.
possible CP-odd Yukawa coupling in the tH process. We fix the sign of the a 1 coupling to be positive, and so the sign of the κ t coupling allows us to float the relative sign of the HVV and Hff contributions in this process. In Scheme 1, we parametrize the tH signal strength with the parameters μ ttH , f Htt CP , and μ V , defined below. There is a weak dependence of the D 0− distributions on f Htt CP in this channel, which we conservatively neglect. We also neglect anomalous HVV couplings in the tH parametrization because their effect is negligible with the current constraints from analysis of the VBF and VH events. In scheme 2, where we do not have a dedicated category for the top quark coupling measurements, we neglect the dependence of the tH process on the HVV couplings because this would have a negligible effect on the results and unnecessarily complicate the fit parametrization.
For EW processes (VBF, ZH, WH), the common unconstrained parameters of interest are μ j ¼ μ V and ⃗ f j ¼ ðf a3 ; f a2 ; f Λ1 ; f Zγ Λ1 Þ in approach 1 or ⃗ f j ¼ ðf a3 ; f a2 ; f Λ1 Þ in approach 2, when using categorization scheme 2 for fitting HVV couplings. We simplify treatment of the EW processes in categorization scheme 1 when fitting for CP-violating Hgg or Htt couplings and allow only the CP-violating f a3 parameter in ⃗ f j . When all anomalous couplings are set to zero, the signal strength μ V is equal to the ratio of the cross sections of all EW processes (VBF, ZH, WH) to the SM expectation. In the case of f ai ¼ 1, this ratio is corrected by the factor ðα ii =α 11 Þ as quoted in Table I for decay cross sections, and the inverse of this factor for the EW production cross sections, due to the evolution of the cross sections with anomalous couplings.
It has been shown that there is no sensitivity to the anomalous couplings in the bbH process [32] and it is parametrized with the signal strength μ j ¼ μ bbH . Depending on the fit implementation discussed in Sec. VI, this signal strength may be correlated with those in other channels, such as μ ggH or μ ttH . The exact treatment of this process is expected to have negligible effect on the results, because it cannot be distinguished kinematically from other dominant processes and its relative contribution in each kinematic region is expected to be negligibly small, as shown in Tables II and III. The background processes i in Eq. (22) include the qq → 4l, gg → 4l, and EW processes, all of which are estimated with simulation, but receive additional constraints from sidebands in data. The EW background includes vector boson scattering and VVV processes, which are the background counterparts of the VBF and VH processes. We also include the ttVV and ttV processes in this background contribution, which are important in the study of the ttH signal process. Interference of the signal and background processes is negligible in the analysis of the on-shell H boson production. The Z þ X background contribution models Z þ jets and other related processes with lepton misidentification and is estimated from the control regions in the data as discussed in Sec. III D. All signal j and background i processes contributing to Eq. (22) and their expected yields are shown in Tables II and III. The signal and background probability distributions P sig jk and P bkg ik appearing in Eq. (22) are binned multidimensional histograms (templates) of observables ⃗ x listed in Table IV. The binning of these templates has been optimized for memory and speed of computer calculations, expected population of events across those bins, and retaining kinematic information. In particular, the large number of discriminants used in the untagged, VBF-2jet, and VH-hadronic categories in scheme 2, requires careful optimization of the binning employed in analysis. In these categories, two bins are used in the two interference discriminants and three bins are used in the other five discriminants, which corresponds to a total of 972 bins in the seven-dimensional (7D) distribution. However, the bins with a very low yield of expected events for all contributions are merged, and the expected symmetry in the distribution of the D CP observable is enforced. As a result, the total number of independent bins depends on the category, but does not exceed 400 in any of the categories. Even though only a limited number of bins is used in each dimension, the 7D distribution retains substantial kinematic information that is nearly optimal for all anomalous couplings targeted in this analysis. This has been validated against a dedicated analysis targeting one anomalous coupling at a time with a much larger number of bins in each dimension for a smaller number of discriminants, as employed in Refs. [16,17]. This nearly optimal performance is realized in large respect due to the optimal population of events across the range of discriminant values by construction of the MELA.
The P sig jk and P bkg ik probabilities depend on the parameters ⃗ ξ jk and ⃗ f j and are therefore interpolated between various templates as a function of these parameters. The ⃗ ξ jk reflect systematic uncertainties either in the normalization or shape of both signal and background templates and an analytical linear interpolation is adopted. The ⃗ f j parameters require nontrivial analytical interpolation of the signal templates, which is discussed in more detail below.

B. Signal parametrization
In the signal production and decay processes, the same anomalous couplings could appear either on both the production and decay sides simultaneously, as in the case of HVV couplings in EW production and H → VV → 4l decay, or only on one side (production or decay). This is illustrated in Eq. (14), with ð P il α ðjÞ il a i a l Þ appearing on the production side and ð P mn α ðfÞ mn a m a n Þ appearing on the decay side. We absorb the width Γ H in Eq. (14) into the overall signal strength and parametrize the kinematic dependence in the signal probability density on the ratio of couplings through ⃗ f j in the following way. In the case of HVV anomalous couplings, we have either L ¼ 5 couplings in approach 1 or L ¼ 4 couplings in approach 2, which we can parametrize with four or three components of ⃗ f j , defined above, and f a1 ¼ ð1 − jf a2 j− jf a3 j − jf Λ1 j − jf Zγ Λ1 jÞ. Let us denote these as f l with l ¼ 1, 2, 3, 4, 5. For the ggH and ttH processes, when we consider anomalous couplings on the production side, we have only two parameters f l ¼ ð1 − jf ggH a3 jÞ and f ggH a3 , or ð1 − jf Htt CP jÞ and f Htt CP , and L ¼ 2. When developing the expression in Eq. (14), one gets a polynomial in the couplings a l ∝ ffiffiffiffiffiffiffi jf l j p signðf l Þ following Eq. (19), which is either quartic (C ¼ 4), in the case of the EW processes, or quadratic (C ¼ 2), when the couplings appear only on the decay or production side. Parametrization of the anomalous coupling dependence of the P sig jk probability density in Eq. (22) is different in these two cases.
This leads to the following general expression for the probability density of the EW processes with C ¼ 4 The following general expression applies to the probability density of the processes with C ¼ 2 In both Eqs. (23) and (24), only the kinematic dependence on ⃗ f j is expressed, while the overall normalization can be absorbed into μ j or accounted for as part of the cross section measurement.
In the general case, there are ðC þ L − 1Þ!=ðC!ðL − 1Þ!Þ terms in either Eqs. (23) or (24). This leads to 70 terms in Eq. (23) when we measure four anomalous HVV couplings in production and decay (L ¼ 5, C ¼ 4), 15 terms in Eq. (24) when we measure four anomalous HVV couplings in decay only (L ¼ 5, C ¼ 2), and three terms in Eq. (24) when we measure Hgg or Htt couplings in production only (L ¼ 2, C ¼ 2). When both HVV (in decay) and Hgg or Htt (in production) couplings are measured at the same time in a given process, the number of terms is multiplied (15 × 3 ¼ 45), since the two probabilities factorize.
The P sig jk;lmðnpÞ templates are extracted from simulation, discussed in Sec. III C, typically using from 3 to 12 samples generated with various ⃗ f j values chosen to map different points of phase-space well and reweighted with the MELA package to cover all possibilities with ðC þ L − 1Þ!= ðC!ðL − 1Þ!Þ combinations of couplings. In parametrizing the signal templates P sig jk;lmðnpÞ , it is important to ensure that the expected number of events in every bin of the probability densities, defined in Eqs. (23) and (24), remains nonnegative at all possible values of ⃗ f j , because a negative yield would cause the likelihood function used for the final fit to become ill defined.
To detect a negative event yield, we first minimize Eqs. (23) or (24), which are polynomials in ffiffiffiffiffiffiffi jf l j p , by finding where the gradient is zero. For Eq. (24), which is quadratic, this is a simple problem in linear algebra. Equation (23) is quartic, so its gradient is a system of cubic equations, which cannot be solved exactly for L > 1. We use the HOM4PS program [105][106][107] to numerically solve the system of equations. If we find the minimum to be negative, we adjust the estimated P sig jk;lmnp until the yield is always positive. This adjustment is made using the statistical uncertainty on P sig jk;lmnp through a cutting planes algorithm [108] implemented using the GUROBI program [109]. In all cases, it is found that only small modifications to the initial estimates of P sig jk;lmnp are needed. The parametrization in Eqs. (23) and (24) is written in the self-consistent full-amplitude approach. In the EFT interpretation of the amplitude fit, the series in powers of jf ai j 1=2 corresponds to terms of different dimension, as discussed in Sec. II C. For example, in Eq. (23) the term with f 2 a1 corresponds to the SM-like contribution with dimension-4 operators, while f 3=2 a1 jf ai j 1=2 corresponds to interference between the SM amplitude and the dimension-6 contributions in the EFT expansion. Assuming that f a1 ∼ 1 and all f ai are small, all other terms could in principle be neglected as a test of EFT validity. In practice, however, neglecting those terms can easily lead to a negative probability in certain points in the phase space of observables, which invalidates the maximum likelihood fit. This does not necessarily mean that the EFT approach is not valid. This happens because the sizable interference terms can lead to a negative partial sum in Eq. (23), especially in the optimized multidimensional space of observables that are sensitive to such interference effects. Therefore, in practice, the fits presented in Sec. VI do not allow one to place constraints without the full series shown in Eqs. (23) and (24).

C. Likelihood fit
The final constraints on parameters μ j and ⃗ f j are placed using the profile likelihood method implemented in the RooFit toolkit [110] within the ROOT [111] framework. The extended likelihood function is constructed using the probability densities in Eq. (22), with each event characterized by the discrete category k and observables ⃗ x. The likelihood L is maximized with respect to the nuisance parameters ⃗ ξ jk describing the systematic uncertainties discussed below, and μ j and ⃗ f j parameters of interest. The allowed 68 and 95% CL intervals are defined using the profile likelihood function, −2Δ ln L ¼ 1.00 and 3.84, for which exact coverage is expected in the asymptotic limit [112].
The reinterpretation of the primary μ j and ⃗ f j results in terms of couplings is performed with the help of Eq. (14) to relate signal strength μ j to couplings and Eqs. (16)-(18) to relate ⃗ f j to coupling ratios. In this way the couplings directly enter the parametrization in Eq. (22). However, without further constraints on the H boson width, such a fit would not provide useful constraints on the coupling size. Therefore, in the total width Γ H parametrization in Eq. (14), we assume that there are no unobserved or undetected H boson decays. We express the width Γ H as the sum of partial decay widths of nine H boson decay modes dominant in the SM, H → bb being the largest and H → μ þ μ − being the smallest. Each partial decay width is scaled as a function of anomalous couplings following the parametrization in Ref. [33].

D. Systematic uncertainties
Several systematic uncertainties are considered in the set of constrained parameters ⃗ ξ jk . The relative expected yields in different categories and the template shapes describing probability distributions in Eq. (22) are varied within either theoretical or experimental uncertainties. All results reported at 68 and 95% CL are dominated by statistical uncertainties. All systematic uncertainties are treated as correlated between different time periods, except for the jetrelated uncertainties, which originate from statistically independent sources, and luminosity uncertainties, which are partially correlated [113][114][115].
The theoretical uncertainties considered are PDF parametrization, factorization, and renormalization scales, the hadronization scale used in PYTHIA, and the underlying event variations. The underlying event modeling uncertainty is determined by varying initial-and final-state radiation scales between 0.25 and 4 times their nominal value. The effects of the modeling of hadronization are determined by simulating additional events with the variation of the nominal PYTHIA tune described in Sec. III. Experimental uncertainties involve jet energy calibration and b-quark-tagging efficiency uncertainties, which are only relevant when production categories are considered, and lepton efficiency and momentum uncertainties, which are similar for the different processes and categories. In the estimation of the Z þ X background, the flavor composition of QCD-evolved jets misidentified as leptons may be different in the Z þ 1l and Z þ 2l control regions, and together with the statistical uncertainty in the Z þ 2l region, this uncertainty accounts for about a AE30% variation in the Z þ X background.
The normalization of the background processes derived from the MC simulation is affected by a 1.8% uncertainty in the integrated luminosity. However, all results are found to be insensitive to luminosity or theoretical constraints on the nonresonant ZZ=Zγ Ã → 4l background. Even though the theoretical cross section is one of the constraints in the fit, the wide sideband included in the range 105 < m 4l < 140 GeV constrains this nearly flat background from the data, with nearly identical results when the theoretical constraints are removed. The main distinguishing feature of this background is the dominant contribution of the Zγ Ã intermediate state, which allows effective separation of the already small background from signal using kinematic information.

VI. RESULTS
The signal strength μ j and the set of parameters ⃗ f j describing the tensor structure of interactions are constrained in each production process and decay H → ZZ=Zγ Ã =γ Ã γ Ã → 4l. In the following, we describe the measurement of f ggH a3 in the ggH process, f Htt CP in the ttH and tH processes, and the combination of the two where the top quark contributes to the gluon fusion loop. We then report measurements of ðf a2 ; f a3 ; f Λ1 ; f Zγ Λ1 Þ in the VBF and VH processes along with the H → 4l decay in all production processes, following approach 1 with the coupling relationship a WW i ¼ a ZZ i . We also report measurements of ðf a2 ; f a3 ; f Λ1 Þ following approach 2 within SMEFT, discussed in Sec. II. These results are interpreted in terms of constraints on the Hff, Hgg, and HVV operators. While all operators could potentially be constrained in a joint analysis of all H boson decay modes, in this paper we analyze only the H → 4l decay mode, and perform a combination with the tH and ttH processes in the H → γγ decay mode. Therefore, for the purpose of illustration, we make further assumptions on how certain couplings, to which this analysis is not sensitive, are related. For example, we must make certain assumptions about the relationship between the Hbb, Hcc, Hττ, and Hμμ couplings and other couplings. These assumptions are discussed in each of the applications presented below.

A. Constraints on Hgg couplings
The measurement of anomalous couplings of the H boson to gluons is presented in Fig. 14 and Table V. Since the direct couplings of the H boson to SM fermions in the gluon fusion loop and to potentially new particles appearing in the loop can not be resolved using this measurement alone, both effects are characterized with two parameters, f ggH a3 and μ ggH . The signal strength μ ggH , which is the ratio A. M. SIRUNYAN et al. PHYS. REV. D 104, 052004 (2021) of the measured cross section of the gluon fusion process to that expected in the SM, is profiled when the f ggH a3 results are reported. The measurement of the f ggH a3 is consistent with zero, as expected in the SM. This can be clearly seen from the D ggH CP and D ggH 0− distributions in Fig. 10. The measured value of μ ggH ¼ 0.86 þ0. 13 −0.11 is consistent with that reported in Ref.
[78] without the fit for the CP structure of interactions. The values of μ ggH and f ggH a3 are uncorrelated. The signal strength of the VBF and VH processes μ V and their CP properties f a3 are also profiled when this measurement is performed. This measurement is also performed simultaneously in a fit with the ttH process with the μ ttH and f Htt CP parameters unconstrained, as discussed below. The tH process is always included with the ttH process with its signal strength expressed through the μ ttH , μ V , and f Htt CP parameters. The parameters f ggH a3 and μ ggH are equivalent to the measurement of the CP-even and CP-odd couplings on the production side, while the HVV couplings on the decay side are constrained from the simultaneous measurement of the VBF and VH processes with f a3 and μ V profiled. The c gg andc gg couplings, introduced in Eqs. (12) and (13), can be extracted from the above measurements. We follow the parametrization o 1f the cross section and the total width from Ref. [33], where c gg ,c gg , κ t ,κ t , κ b , andκ b contribute. Since it is not possible to disentangle all these couplings in a single process, we fix κ t  SM expectation and leave c gg andc gg , which describe possible BSM contribution in the loop, unconstrained. The small contribution of the H → γγ and Zγ decays to the total width is assumed to be SM-like. The resulting constraints on c gg andc gg are shown in Fig. 14, right. The general features of these constraints are the following. The pure signal strength measurement μ ggH , available even without the fit for f ggH a3 , provides a constraint in the form of a ring on a two-parameter plane in Fig. 14, right. The measurement of f ggH a3 resolves the areas within this ring. Since the sensitivity of the f ggH a3 measurement is currently just under 68% CL, this resolution is not strong. The H boson width dependence on c gg andc gg is relatively weak and does not alter this logic considerably. The results are consistent with the SM expectation.
As mentioned earlier, it is not possible to resolve the loop contributions from the SM or BSM particles in this measurement. Therefore, the deviations of the SM-like Yukawa couplings κ t and κ b from unit values are absorbed into the effective c gg measurement, and the CP-odd Yukawa couplingsκ t andκ b are absorbed into the effectivec gg measurement, together with possible contributions from BSM particles. However, reinterpretation of these results is possible in terms of the independent Yukawa couplings and effective pointlike gluon couplings in combination with the ttH and tH modes, as discussed below.

B. Constraints on Htt couplings
The measurement of anomalous couplings of the H boson to top quarks is presented in Fig. 15 and Table V. First, the measurements of f Htt CP from the ttH and tH processes only are reported. The signal strength μ ttH , which is the ratio of the measured cross section of the ttH process to that expected in the SM, is profiled when the f Htt CP results are reported. The measured value of μ ttH ¼ 0.17 þ0. 70 −0.17 is consistent with that reported in Ref.
[78] without the fit for the CP structure of interactions. In both cases we observe downward fluctuations in the signal yield compared to expectation, but these fluctuations are not statistically significant. There is no significant linear correlation between μ ttH and f Htt CP . The signal strength of the VBF and VH processes μ V , ggH process μ ggH , and their CP properties f a3 and f ggH a3 are also profiled when this measurement is performed. This analysis of the ttH and tH processes is not sensitive to the sign of f Htt CP . However, for later combination with the ggH measurement, presented above, under the assumption of the top quark dominance in the gluon fusion loop, symmetric constraints on f Htt CP are reported.
The observed best fit f Htt CP value gives preference to the CP-odd Yukawa coupling. This comes from the negative value of the D ttH 0− discriminant for the one observed signallike event in Fig. 9. However, this result is statistically consistent with the pure CP-even Yukawa coupling expected in the SM. With just about two signal ttH events and H → 4l decays. Left: observed (solid) and expected (dashed) likelihood scans of f Htt CP in the ttH and tH processes in the H → 4l (red), γγ (black), and combined (blue) channels, where the combination is done without relating the signal strengths in the two processes. The dashed horizontal lines show 68 and 95% CL Right: observed confidence level intervals on the κ t andκ t couplings reinterpreted from the f Htt CP , μ ttH , and μ V measurements in the combined fit of the H → 4l and γγ channels, with the signal strengths in the two channels related through the couplings as discussed in text. The dashed and solid lines show the 68 and 95% CL exclusion regions in two dimensions, respectively. and many fewer tH events expected to appear in the fit in the H → 4l channel under the assumption of the SM cross section, according to Table II, the expected confidence sensitivity on the f Htt CP constraints is low. Nonetheless, the high signal purity in the H → 4l channel implies that every observed event candidate carries a large statistical weight. The importance of including the CP measurements in the ttH and tH production modes also becomes evident when combination with ggH is performed. There is a significant gain in such a combination beyond a simple addition of independent measurements, as discussed in Sec. VI C.
The CMS experiment recently reported the measurement of the f Htt CP parameter in the ttH and tH production processes with the decay H → γγ [26] (shown also in Table V). In that measurement, the signal strength μ γγ ttH parameter is profiled, while the signal strengths in other production processes are fixed to the SM expectation. However, there is a very weak correlation of the measurement in the ttH and tH processes with parameters in the other production mechanisms. Therefore, we proceed with a combination of the f Htt CP measurements in the H → 4l and γγ channels, where we correlate their common systematic uncertainties, but not the signal strengths of the processes. In particular, we do not relate the μ ttH and μ γγ ttH signal strengths because they could be affected differently by the particles appearing in the loops responsible for the H → γγ decay. The results of this combination are presented in Fig. 15 and Table V. The measured signal strength in the H → 4l channel is μ ttH ¼ 0.04 þ0.76 −0.04 , uncorrelated with f Htt CP , while μ γγ ttH ¼ 1.23 þ0.33 −0.24 and the correlation with f Htt CP is þ0.20. The pure pseudoscalar hypothesis of the H boson corresponding to f Htt CP ¼ 1 in the case of the CP-odd Yukawa interaction is excluded at 3.2 standard deviations, while the expected exclusion is 2.7 standard deviations. Below, we also present an interpretation of these results where the signal strengths in the two H boson decay channels are related through the couplings.
In the above measurements, the f Htt CP parameter has the same meaning in both the H → 4l and γγ channels. In order to make an EFT coupling interpretation of the results, we have to make a further assumption that no BSM particles contribute to the loop in the H → γγ decay. Without this or a similar assumption, the signal strength in the H → γγ decay cannot be interpreted without ambiguity. We further reparametrize the cross section following Ref.
[33] with the couplings κ t andκ t , and fix κ b ¼ 1 and κ b ¼ 0. The bottom quark coupling makes a very small contribution to the loop in the H → γγ decay, but it makes a large contribution to the total decay width, where we assume that there are no unobserved or undetected H boson decays. In order to simplify the fit, we do not allow anomalous HVV couplings, and the measurement of the signal strength μ V constrains the contribution of the a 1 coupling in the loop. The f ggH a3 and μ ggH parameters are profiled in this fit.
The observed confidence level intervals on the κ t andκ t couplings from the combined fit of the H → 4l and γγ channels are shown in Fig. 15. There is no linear correlation between the values of κ t andκ t . As was the case for the (c gg ,c gg ) measurement in Fig. 14, the pure yield measurement in the ttH process would constrain a ring in the two-dimensional plane. However, the CP-sensitive measurement of f Htt CP disfavors the values away fromκ t ¼ 0. Moreover, the sign ambiguity between the κ t and −κ t values cannot be resolved in the ttH channel alone. With the inclusion of the tH process, the negative values of κ t are disfavored because strong constructive interference between the amplitudes induced by the HVV and Htt couplings would result in enhanced tH yield, inconsistent with the data. Therefore, the sign of κ t is defined in reference to the tree-level HVV coupling a 1 . But, the sign ambiguity between theκ t and −κ t values cannot be resolved in this fit, unless information from the other channels is incorporated, such as information from the gluon fusion loop discussed below.

C. Constraints on Htt and Hgg couplings in combination
First, we consider the ggH process under the assumption of top quark dominance in the gluon fusion loop. The measurement of anomalous couplings of the H boson to top quarks for this case is presented in Fig. 16 and Table V. Similar to the case of the H → γγ loop discussed above, the cross section of the ggH process, normalized to the SM expectation, is parametrized following Ref.
[33] to account for CP-odd Yukawa couplings as follows: where we set κ f ¼ κ t ¼ κ b andκ f ¼κ t ¼κ b . Equation (25) sets the relationship between f Htt CP and f ggH a3 , reported in Fig. 14 and Table V, according to Eq. (17).
Constraints on f Htt CP are also shown for the combination of the ttH, tH, and ggH processes with H → 4l in Fig. 16 and Table V. The combination of the H → 4l and γγ channels with ttH, tH, and ggH processes proceeds in a similar manner and is also shown in Fig. 16 and Table V. In this case, we do not allow anomalous HVV couplings, and the measurement of the signal strength μ V constrains the contribution of the a 1 coupling in the H → γγ loop.
The gain in this combination of the ggH and tH and ttH processes is beyond the simple addition of the two constraints. While in the ggH and ttH analyses the signal strengths of the two processes are independent, they could be related under the assumption of top quark dominance in the loop using Eq. (25). As discussed in Sec. II, CP-odd coupling predicts rather different cross sections in the two processes: σðκ f ¼ 1Þ=σðκ f ¼ 1Þ is 2.38 in the gluon fusion process dominated by the top quark loop and 0.391 in the ttH process. This means that the ratio differs by a factor of 6.09 for f Htt CP ¼ 1 when compared to SM (f Htt CP ¼ 0). This correlation enhances the sensitivity in the f Htt CP measurement. For example, the combined sensitivity from tH and ttH (with either H → 4l alone or together with H → γγ), and ggH is significantly improved compared to separate analyses, and the result is not just a simple addition of two independent results, as shown in Fig. 16 and Table V. This effect also enhances correlation between f Htt CP and the yield parameters. In the full combination, the measured signal strength is μ ttH ¼ 0.70 þ0.30 −0.25 and the correlation with f Htt CP is þ0.96. Finally, we present the reinterpretation of the f ggH a3 , f Htt CP , and signal strength measurements in terms of constraints on c gg ,c gg , κ t , andκ t . In this fit, it is assumed that κ b ¼ κ c ¼ κ μ ¼ 1 andκ b ¼κ c ¼κ μ ¼ 0 in the fermion coupling contribution to the loops and in the decay width parametrization [33]. The gluon fusion loop is parametrized in terms of pointlike couplings c gg andc gg , and the top and bottom quark contributions. These pointlike and top couplings can be resolved in combination with the ttH and tH production processes. The H → γγ loop is parametrized with the top and bottom quark and with the W contributions. One cannot generally relate the pointlike couplings in this loop and in the gluon fusion loop, and they are assumed to be zero in H → γγ. The measurement of the signal strength μ V constrains the contributions of the a 1 coupling, affecting the W contribution to the H → γγ loop and tH process, and anomalous HVV couplings are not allowed. By convention, the a 1 coupling is constrained to be positive, which sets the relative sign of κ t . It is assumed that there are no unobserved or undetected H boson decays.
The constraints are shown in Fig. 17, where the likelihood scans are plotted for the pairs of parameters: (c gg ,c gg ), (κ t ,κ t ), and (κ t , c gg ), with the other two parameters either profiled or fixed to the SM expectation. On the likelihood scan of (c gg ,c gg ) with the other parameters fixed, the appearance is generally similar to the scan shown in Fig. 14. This is because the addition of the ttH and tH processes does not alter the results much with their couplings fixed to the SM. However, with the other parameters affecting the gluon fusion loop left floating, the contours are washed out, as one would expect with more degrees of freedom in a fit.
On the likelihood scan of (κ t ,κ t ) with the other parameters profiled, one can observe the effects similar to that in Fig. 15. This is because the ggH process does not bring additional constraints due to uncertainties with the pointlike interactions, but may introduce additional modifications to the fit parameters. However, with the pointlike interactions c gg andc gg set to zero, constraints on the κ t and κ t couplings appear tighter as a result of the combination of information from the ttH, tH, and ggH processes. In both cases, the general features remain similar to Fig. 15, such as pure yield measurements leading to constraints within a ring, CP-sensitive measurements resolving areas within this ring, and the tH process leading to the exclusion of negative values of κ t . However, in this case, ambiguity between the positive and negative values ofκ t can be resolved with the inclusion of the ggH process, where the D ggH CP discriminant carries information sensitive to the sign. On the likelihood scan of (κ t , c gg ) with the other parameters fixed, we observe a resolved fourfold ambiguity of the best fit ranges. Within each range, there is a large correlation between the two parameters. This happens because the pointlike interaction c gg is equivalent to a BSM heavy quark Q contribution to the loop. It is hard to distinguish between such a new heavy quark and the heavy top quark. The two amplitudes add constructively, leading to a large anticorrelation. The rate of the gluon fusion process would be roughly proportional to ðκ t þ ακ Q Þ 2 with α ∼ 1. Given that the ttH rate constrains κ t ∼ AE1, the ggH rate would constrain (κ t , κ Q ) to four discrete sets of values around ðþ1; 0Þ, ðþ1; −2Þ, ð−1; 0Þ, and ð−1; þ2Þ. The presence of the tH process shifts the preferred negative  (dashed) likelihood scans of f Htt CP are shown in the ggH process with H → 4l (black), ttH, tH, and ggH processes combined with H → 4l (red), and in the ttH, tH, and ggH processes with H → 4l and the ttH and tH processes with γγ combined (blue). Combination is done by relating the signal strengths in the three processes through the couplings in the loops in both production and decay, as discussed in the text. The dashed horizontal lines show 68 and 95% CL exclusion. κ t solutions away from −1 and makes it less likely than the þ1 value for the reasons discussed above. However, the local minimum near κ Q ∼ −2, corresponding to c gg ∼ −0.017, cannot be excluded, even though the global minimum is at c gg ¼ −0.001, close to the null SM expectation. In the case with the other parameters profiled, the constraints on the (κ t , c gg ) plane get washed out further, as expected in a fit with more degrees of freedom. In this case, the CP-odd amplitudes can compensate for some effects of the CP-even ones. However, some sensitivity is retained because CP-sensitive measurements constrain the relative contribution of CP-odd amplitudes.

D. Constraints on HVV couplings
The measurement of anomalous couplings of the H boson to EW vector bosons in approach 1 with the relationship a WW i ¼ a ZZ i is presented in Fig. 18 and Table VI. Figure 18 shows the observed and expected likelihood scans in the simultaneous measurement of f a3 , f a2 , f Λ1 , and f Zγ Λ1 , where the CP-sensitive parameter f ggH a3 and the signal strength parameters μ V and μ ggH are profiled, and where we relate μ ttH and f Htt CP to μ ggH and f ggH a3  assuming top quark dominance in the loop. The results are shown for each coupling separately, with the other three anomalous couplings either set to zero or left unconstrained in the fit. Figure 19 shows the same results presented as two-dimensional contours, where all couplings discussed above are left unconstrained. In all cases, the likelihood scans are limited to the physical range of jf a3 j þ jf a2 j þ jf Λ1 j þ jf Zγ Λ1 j < 1. There are several features visible on these plots. First, the results with all other couplings constrained to zero exhibit narrow minima near f ai ¼ 0 in both the expected and the observed scans. This effect comes from utilizing production information. The anomalous coupling terms in Eq. (2) are multiplied by a factor of q 2 i , which is larger in VBF and VH production than in H → 4l decay. As a result, the cross section in VBF and VH production increases quickly with f ai . At the same time, the constraints above f ai ∼ 0.02 are dominated by the decay information from H → 4l. However, when all four anomalous couplings are allowed to float independently, the best fit value is ðf a3 ; f a2 ; f Λ1 ; f Zγ Λ1 Þ ¼ ð−0.00805; −0.24679; 0.18629; −0.02884Þ. This global minimum is driven by the decay information from H → 4l and is only slightly preferred to the local minimum at (0, 0, 0, 0), with a difference in −2 lnðLÞ of 0.05 between the SM value and the global minimum. The local minimum at (0, 0, 0, 0) is still evident in the four-dimensional distribution and its projections on each parameter, and is driven by the production information, as discussed above for TABLE VI. Summary of constraints on the anomalous HVV coupling parameters with the best fit values and allowed 68 and 95% CL intervals. Three scenarios are shown for each parameter: with three other anomalous HVV couplings set to zero (first), with three other anomalous HVV couplings left unconstrained (second), in approach 1 with the relationship a WW i ¼ a ZZ i in both cases; and with two other anomalous HVV couplings left unconstrained (third), in approach 2 within SMEFT with the symmetry relationship of couplings set in Eqs. (3)- (7). The f Zγ Λ1 parameter is not independent in the latter scenario.
Parameter Scenario Observed Expected f a3 8 > > > > > > > > > > > < > > > > > > > > > > > : Owing to what appears to be a statistical fluctuation in the observed data when the −2 lnðLÞ minima obtained from the decay and from the production kinematics differ, the observed constraints appear weaker than expected. However, the results are still statistically consistent with the SM and with the expected constraints in the SM. Should the global minimum nonetheless persist away from (0, 0, 0, 0) with more data, it will be interesting to study consistency of the constraints from the VBF and VH production and from the H → 4l decay. The production and decay test different ranges of q 2 i , as discussed above. If the q 2 i growth is truncated in the VBF and VH production due to lower-energy BSM effects, then the decay information becomes more important. E. Constraints on HVV couplings within SUð2Þ × Uð1Þ symmetry The above studies of the H boson couplings to EW vector bosons are repeated following approach 2 within SMEFT with the symmetries in Eqs. (3)- (7). In this case, the f Zγ Λ1 parameter is not independent. Therefore, constraints on the three parameters f a3 , f a2 , f Λ1 , and the signal strength are obtained in this scenario following the same approach as above. These constraints are shown in Fig. 20 and Table VI. The measured signal strength is TABLE VII. The observed correlation coefficients of the signal strength μ V and the f a3 , f a2 , f Λ1 parameters in approach 2 within SMEFT with the symmetry relationship of couplings set in Eqs. (3)-(7).

Parameter
Observed correlation  Table VII. Keeping only linear terms and dropping terms with order greater than one for anomalous couplings does not allow us to make a reasonable likelihood scan, since the probability density goes negative, as discussed in Sec. II C.
Since the relationship of the HWW and HZZ couplings does not affect the measurement of the f a3 parameter in the H → 4l decay, the constraints from the decay information in the wider range of f a3 in approach 2 are unaffected compared to approach 1, when other couplings are fixed to zero. However, with one less parameter to float, the constraints are modified somewhat when all other couplings are left unconstrained. The modified relationship between the HWW and HZZ couplings also leads to some modification of constraints using production information in the narrow range of f a3 . On the other hand, the f a2 and f Λ1 parameters are modified substantially because the f Zγ Λ1 information gets absorbed into these measurements through symmetry relationships.
The measurement of the signal strength μ V and the f a3 , f a2 , f Λ1 parameters can be reinterpreted in terms of the δc z , c zz , c z□ , andc zz coupling strength parameters. Observed one-and two-dimensional constraints from a simultaneous fit of SMEFT parameters are shown in Figs. 21 and 22. The c gg andc gg couplings are left unconstrained. A summary of all constraints on the Htt, Hgg, and HVV coupling parameters in the Higgs basis of SMEFT, including the correlation coefficients, is shown in Table VIII. The results in this table are taken from Secs. VI C and VI E, as measured in the tH, ttH, ggH, and EW processes.
The above interpretation of HVV results in terms of the δc z , c zz , c z□ , andc zz couplings can be extended into an interpretation in terms of the couplings in the Warsaw SMEFT basis [28]. In this basis, nine operators are considered: c H□ , c HD , c HW , c HWB , c HB , c HW , c HWB , c HB , and δ v , where the latter is a linear combination of additional Warsaw basis operators [71]. However, not all nine of these operators are independent. First of all, consideration of Eq. (3) leads to δ v expression as a linear combination of c HD and c HWB . Four constraints on the couplings a γγ;Zγ 2 and a γγ;Zγ 3 lead to only one of the three operators c HW , c HWB , and c HB being independent, and only one of c HW , c HWB , and c HB being independent. Therefore, we obtain only four independent constraints, the same number as in the Higgs basis. We note that the couplings of the Z boson to fermions are fixed to those expected in the SM because those are well constrained from prior measurements and this constraint is already included in our primary measurements. Even though some of the above EFT operators may affect couplings of the Z boson, their effect must be compensated by the other EFT operators not affecting the H boson couplings directly. With the above constraints, we use the tools in Refs. [33,71] to relate TABLE VIII. Summary of constraints on the Htt, Hgg, and HVV coupling parameters in the Higgs basis of SMEFT. The observed correlation coefficients are presented for the Htt and Hgg and HVV couplings in the fit configurations discussed in text and shown in Figs operators in the Higgs and Warsaw bases. Since it is arbitrary which one of the three operators is considered to be independent, we present results with all three choices in each case. The results can be found in Table IX, where three other independent couplings are left unconstrained for each measurement.

VII. SUMMARY
In this paper, a comprehensive study of CP violation, anomalous couplings, and the tensor structure of H boson interactions with electroweak gauge bosons, gluons, and fermions, using all accessible production mechanisms and the H → 4l decay mode, is presented. The results are based on the 2016-2018 data from pp collisions recorded with the CMS detector during run 2 of the LHC, corresponding to an integrated luminosity of 137 fb −1 at a center-of-mass energy of 13 TeV. These results significantly surpass our results from run 1 [13] in both precision and coverage. The improvements result not only from a significantly increased sample of H bosons but also from a detailed analysis of kinematic distributions of the particles associated with the H boson production in addition to kinematic distributions in its decay. These results also surpass our earlier studies of on shell production of the H boson in this decay channel with a partial run 2 dataset [16,17].
The parametrization of the H boson production and decay processes is based on a scattering amplitude or, equivalently, an effective field theory Lagrangian, with operators up to dimension six. Additional symmetries and prior measurements allow us to reduce the number of independent parameters and make a connection to the standard model effective field theory (SMEFT) formulation. Dedicated Monte Carlo programs and matrixelement reweighting techniques provide modeling of all kinematic effects in the production and decay of the H boson, with any variation of parameters of the scattering amplitude and with full simulation of detector effects. Each production process of the H boson is identified using the kinematic features of its associated particles. The MELA is employed to construct observables that are optimal for the measurement of the targeted anomalous couplings in each process, including CP-sensitive observables. A maximum likelihood fit allows a simultaneous measurement of up to five HVV, two Hgg, and two Htt couplings.
For the first time, we present a complete and dedicated study of CP properties in the H boson coupling to gluons through a loop of heavy particles using the CP-sensitive observables, while separating the electroweak and strong boson fusion processes. An interpretation of the loop contribution is made both with and without an assumption of top quark dominance, which allows for a new heavy particle to contribute. In both cases, combination with the CP-sensitive measurement of the Htt coupling in the ttH and tH processes allows either simultaneous or separate measurements of the two effective pointlike Hgg couplings and the two Htt couplings, both CP-odd and CP-even. For the ttH and tH processes, results in the H → 4l channel are combined with those in the H → γγ channel [26]. This is the first comprehensive study of CP properties in the Hgg and Htt couplings from a simultaneous measurement of the ggH, ttH, and tH processes.
Also for the first time, we present the measurement of CP properties and the tensor structure of the H boson's interactions with two electroweak bosons with up to five parameters measured simultaneously. The HVV coupling is analyzed in VBF and VH production and in H → VV → 4l decay. The measurements are performed with two approaches. In the first approach, more lenient symmetry considerations are applied, which allow a less restrictive interpretation of results. In the second approach, SUð2Þ × Uð1Þ symmetry is invoked and the formulation becomes equivalent to SMEFT. The operator basis is chosen to coincide with the couplings of the mass eigenstates, which allows us to minimize the number of independent parameters. A translation of the SMEFT results to the Warsaw basis is also presented for easier comparison with other results.
In all cases, we first present results in terms of the total cross section of a process and the fractional contribution of each anomalous coupling. These results are further reinterpreted in terms of direct constraints on the couplings by applying certain assumptions about the H boson total width. Each of the measurements presented here is limited by statistical precision and is consistent with the expectations for the standard model Higgs boson.  [35] ATLAS Collaboration, Measurements of the Higgs boson inclusive and differential fiducial cross sections in the 4l decay channel at ffiffi ffi s p ¼ 13 TeV, Eur. Phys. J. C 80, 942 (2020).