Constraining anomalous Higgs boson couplings to the heavy flavor fermions using matrix element techniques

In this paper we investigate anomalous interactions of the Higgs boson with heavy fermions, employing shapes of kinematic distributions. We study the processes $pp \to t\bar{t} + H$, $b\bar{b} + H$, $tq+H$, and $pp \to H\to\tau^+\tau^-$, and present applications of event generation, re-weighting techniques for fast simulation of anomalous couplings, as well as matrix element techniques for optimal sensitivity. We extend the MELA technique, which proved to be a powerful matrix element tool for Higgs boson discovery and characterization during Run I of the LHC, and implement all analysis tools in the JHU generator framework. A next-to-leading order QCD description of the $pp \to t\bar{t} + H$ process allows us to investigate the performance of MELA in the presence of extra radiation. Finally, projections for LHC measurements through the end of Run III are presented.


I. INTRODUCTION
The discovery [1,2] of the H boson by the ATLAS and CMS experiments during Run I of the LHC marked an important milestone in the evolution of our understanding of fundamental particle physics. One of the most important goals now is a precise understanding of the newly discovered state, including its couplings to other particles as well as its CP nature. Any significant deviation from Standard Model (SM) predictions would reveal the existence of new physics in the Higgs sector and should to be classified according to its anomalous coupling structures. Likewise, in the case of a discovery of a new resonance at the LHC, a similar program of investigating properties and couplings is required.
Since experimental efforts during Run I mostly focused on the H boson decaying to a pair of vector bosons 1 , extensive studies of the HV V couplings and corresponding CP properties were performed [3][4][5][6], leading to results consistent with the Standard Model nature of the H boson with quantum numbers J P C = 0 ++ . However, many generic models of new physics predict deviations which are beyond the current experimental precision [7]. This leaves ample room for anomalous interactions to hide as small modifications of the SM structure. Moreover, a complete knowledge of the SM Higgs mechanism requires the study of the H boson interactions with fermions. While the observation of Hgg and Hγγ interactions established the Htt coupling through a closed loop, a detailed understanding only arises from the observation of various mass hierarchies, as well as the minimal flavor-universal Yukawa interaction as predicted by the SM. Hence, it is of paramount importance to investigate possible anomalous coupling patterns with which different quarks and leptons may interact with the Higgs field. The most promising approach is the study of differential distributions in processes with direct sensitivity through the associated production of the H boson with on-shell fermions, f f + H, or through the decay H → ff .
There has been considerable effort in modeling the Hff couplings and developing tools for their analysis in associated production  and in the H → τ + τ − decay [32][33][34][35][36][37][38]. Current experimental analyses have measured the coupling strength of Hbb and Htt only through closed loops [4,6]. There have been experimental searches for the H boson production in association with a single top quark [39] and with tt [40][41][42][43], with strong evidence for the latter in Run I of the LHC. The process bbH has not been studied with a dedicated analysis so far but there is evidence for the H boson decay into bb pairs [44][45][46]. The decay of H → τ + τ − is observed when results of CMS [47] and ATLAS [48]

II. PARAMETERIZATION OF HIGGS BOSON COUPLINGS
We describe the interactions between a spin-zero particle H and two fermions through the amplitude whereψ f and ψ f are the Dirac spinors, m f is the fermion mass, and v is the SM Higgs field vacuum expectation value. In the SM, the couplings 3 have the values κ f = 1 andκ f = 0. Any deviation from these values indicates the presence of physics beyond the SM, which may for example arise through heavy loop-induced fields. In particular, theκ f coupling parameterizes the contribution of a CP -odd pseudoscalar boson, and CP violation occurs when both κ f andκ f are nonzero. One may equivalently choose to express the couplings through a Lagrangian (up to an unphysical global phase) which allows a connection to be made between the couplings κ f andκ f and anomalous operators in an effective field theory. We assume the couplings to be independent of kinematics, which corresponds to dimension-six operators in the effective field theory. Higher-dimension operators could easily be considered through q 2 -dependent couplings in our framework, where q is the momentum transfer. However, in our study we neglect these higher-dimension contributions because they are expected to be small. The hermiticity of the Lagrangian requires κ f andκ f to be real. Nevertheless, in order to consider the broadest range of scenarios, we allow κ f andκ f to be complex, and trust that, should the unitarity of scattering amplitudes be violated as a result, it will be restored in the full theory. It is convenient to parameterize anomalous couplings through a mixing angle, with κ f ∝ cos α andκ f ∝ sin α. Equivalently, we introduce the fractions where the f CP parameter is conveniently bounded between 0 and 1, is uniquely defined, and can be interpreted as the cross section fraction corresponding to the pseudoscalar coupling, and therefore is directly related to experimentally observable effects. It is a convenient counterpart of the f a3 parameter defined for the HV V couplings [3,5,63]. While the phase φ CP can in general take any value between 0 and 2π, it is reasonable to assume that the ratioκ f /κ f is real, that is φ CP = 0 or π. However, we do not need to impose this restriction and will also consider other values of φ CP . The parameters f CP and φ CP in principle depend on the fermion couplings under consideration and should be denoted f Hff CP and φ Hff CP , but in most cases this will be clear from the context. The tqH production also involves the HW W coupling. We therefore recall the coupling of the H to two vector bosons [63] where i is the polarization of the vector boson of mass m V and momentum q i , the field strength tensor is

III. MONTE CARLO SIMULATION
The JHU generator framework [61][62][63] involves both the Monte Carlo generation of unweighted events and the MELA package used in the analysis of the H boson couplings. For top quark pair production in association with a spin-zero boson H we compute the leading-order processes gg → tt + H and qq → tt + H, followed by spin-correlated top quark decays t → bW (→ f f ) in the narrow-width approximation. Any leptonic or hadronic decay mode of the top quarks can be described. Representative Feynman diagrams are shown in Fig. 1, where we allow for the anomalous Htt couplings shown in Eq. (1). The H boson is considered stable in the respective matrix elements describing production, and its decay into any possible channel can be introduced subsequently through processing the generated events using the JHU generator framework.
Since hadronic production of tt + H final states involves color flow in initial and final states, additional jet radiation plays an important role in the description of this process. In fact, almost 40 % of all tt events are accompanied by jets with transverse momentum harder than 40 GeV [66]. It is therefore important to study the impact of radiative corrections on event kinematics and the matrix element observables. To this end, we also calculate the next-to-leading order QCD correction to the pp → tt + H process.
The framework for NLO QCD computations is an extension of the TOPAZ code which two of us developed for anomalous coupling studies of tt + Z final states [67,68]. We calculate the virtual correction to the gg and qq initial states using the numerical implementation of D-dimensional generalized unitarity techniques [69][70][71][72]. The real emission corrections involve the partonic processes gg → ttg, qq → ttg, qg → ttq andqg → ttq, which we regularize using the massive dipole subtraction techniques of Refs. [73,74]. A consistent expansion in the strong coupling constant also requires the computation of the NLO corrections to the top quark decay t → bW and the subsequent W → jj decay. We account for these contributions in the narrow-width approximation using the implementation developed in Ref. [75]. Non-resonant and off-shell effects are expected to scale parametrically as Γ t /m t ≈ 1 % and hence can be safely neglected provided that phase space cuts do not severely constrain the top quark invariant mass. This has been explicitly confirmed in studies for tt + H production at LO [76] and NLO QCD [77].
We obtain the pp → bb + H process from pp → tt + H by replacing m t → m b in the matrix elements and phase space (and removing the top quark decay), while preserving the five flavor scheme with massless initial state quarks. Hence, we neglect the newly appearing t-channel diagram in bb → bb reactions which is, however, doubly-suppressed by the small b quark parton luminosity. In this way the H boson is always emitted from the massive final state quarks only. We believe these approximations are sufficient for studying anomalous interactions in our analysis. Simulation of the single top quark production process in association with a spin-zero boson relies on the partonic processes qb → q t + H (t-channel process) and qq → tb + H (s-channel process). The former topology is shown in Fig. 2 (a) and (b), and the latter topology is shown in Fig. 2 (c) and (d). We make use of analytic expressions for the leading-order SM matrix elements [78] and extend them to include anomalous couplings, keeping the five flavor scheme so that the H boson is never radiated off the initial state b-quark. An extensive comparison of the four and five flavor schemes for this process was performed in Ref. [26]. Here, we only note that differences between the two schemes are due to missing higher orders in the truncated perturbative series. Since the perturbative convergence for this process is good, these ambiguities are adequately captured by the scale uncertainty. Interestingly, in contrast to ttH and bbH production, tqH production includes not only Hff coupling but also HW W coupling (depicted in Fig. 2 (b) and (d)). The interference between these diagrams is destructive and leads to a strongly suppressed production rate in the SM [79]. Therefore any new physics modification of either the Hff or HW W coupling may spoil this suppression and lead to a substantially enhanced production rate and altered kinematics. We therefore include anomalous HW W couplings, following Eq. (4).
The study of spin-zero H boson anomalous couplings to tau leptons relies on the matrix element H → τ + τ − with subsequent spin-correlated decays τ → µ ν τ ν µ or τ → q q ν τ . These decay chains are illustrated in Fig. 3. This decay mode supplements the existing H → V V decays (V = W, Z, γ, g) within the JHU generator framework such that any H boson production process can be interfaced with this decay. An option for stable τ leptons allows one to study H → µ + µ − or H → e + e − decays as well. Moreover, the tau decay chain encompasses the same structure as the top quark decay, enabling the future study of fully spin-correlated decay pp → X → tt → bbW (→ f f )W (→ f f ), where X is any massive spin-zero state. In this work, we will only focus on anomalous coupling studies in pp → H → τ + τ − . In the current implementation, the form factors for hadronic tau decay are not implemented in the generator and instead the inclusive tau decay is simulated. Below we illustrate the re-weighting technique to obtain the H → τ + τ − process with anomalous couplings using SM simulation of hadronic τ decays with hadronic form factors by the TAUOLA package [80].
The generation of unweighted events for H boson production in association with heavy-flavor quarks is performed at leading order (LO), complemented with parton shower generated by PYTHIA8 [81,82]. The H boson decay is simulated independently from its production. In all cases, the Les Houches Event (LHE) file format [83] is used to interface the JHU generator program. We also generate weighted events at NLO in QCD for ttH production to investigate the impact of radiative corrections, as discussed above. The simulation of the SM processes ttH and bbH has been checked against the NLO QCD production simulation by POWHEG [84][85][86], pseudoscalar ttH production at NLO QCD has been checked against Ref. [87], while the H → τ + τ − decay is validated with TAUOLA. The background ttV V samples in this study are generated with MadGraph [88].
In the following, we focus on the LHC energy of   as well as the NNPDF3.0 parton distribution functions [89]. We summarize the processes relevant to our study of H boson CP properties in heavy flavor fermion interactions, discussed above, in Table I. We also show their SM production cross sections and the order in perturbative QCD to which they are simulated. For each process shown, we provide the matrix elements through the MELA library. One direct application of MELA is kinematic discriminants for the optimal analysis, as discussed in Section IV. This technique also allows one to re-weight an existing Monte Carlo sample to any model with anomalous couplings without the need for additional CPU-consuming simulation. This is particularly important for the LHC experiments where modeling ATLAS and CMS detector response sometimes requires months of wall-clock time. A successful application of this procedure has been presented in Ref. [3].

IV. MATRIX ELEMENT TECHNIQUE
The matrix elements, or multivariate per-event likelihoods, maximize the amount of information that can be extracted from a given event. These techniques were used for example in top and bottom quark, as well as electroweak boson measurements, and proved to be powerful tools for the H boson discovery and characterization during Run I of the LHC on both CMS [2,3,41,52,53,55,57,91] and ATLAS [5,59,60] experiments. As part of the latter development, we investigated application of these techniques to the production and decay processes involving H boson coupling to vector bosons in Refs. [61][62][63]. Here we extend the MELA technique to the processes involving H boson coupling to heavy flavor fermions.
We take the gg(qq) → ttH processes as an example to define a complete set of kinematic observables following the full sequence of the process, similar to the HV V production and decay kinematics discussed in Refs. [61,63]. These observables are equivalent to the more familiar observables defined in the laboratory frame, as shown in Appendix A, but provide a more intuitive insight into the production and decay dynamics. We then define the complete set of matrix element discriminants, following Refs. [62,63], in application to the processes involving heavy fermion couplings to the H boson.

A. Kinematics in the H boson production and decay
The processes gg(qq) → ttH, tqH, or bbH with subsequent decay of the top quarks and the H boson can be characterized by the four-momenta of the decay products, such as leptons and quark jets. In the case of one final state neutrino, its momentum can be deduced from a kinematic fit using mass constraints and utilizing the missing transverse energy information. In the following description we consider the ttH production in its center-of-mass frame. Both longitudinal and transverse momenta of the ttH system can be parameterized separately. They are driven by QCD effects, either parton distribution functions of the proton for rapidity or additional jet radiation for transverse momentum.
Similar to the description of the H boson production and decay with couplings to vector bosons [61][62][63], it is convenient to describe the complete kinematics of the process by a set of angles and invariant masses, which we generically denote as Ω, following the sequential processes. The definition of observables in the process gg(qq) Fig. 4. The following set of angles and invariant masses is  • θ H : angle between the H boson direction and the incoming partons in the ttH frame; • θ * V : angle of the H → V V (ff ) decay with respect to the opposite tt direction in the H frame; • Φ * V : angle between the production plane, defined by incoming partons and H, and H → V V (ff ) decay plane; • θ t : angle between the top quark direction and the opposite Higgs direction in the tt frame; • Φ * t : angle between the decay planes of the tt system and H → V V (ff ) in the ttH frame; • m tt : invariant mass of the tt system; • θ W : angle between W + and opposite of the bb system in the W + W − frame; • Φ W : angle between the production (bb)(W + W − )H plane and the plane of the W + W − system in the tt frame; • θ b : angle between the b quark and opposite of the W + W − system in the bb frame; • Φ b : angle between the planes of the bb and W + W − systems in the tt frame; • m W b1 or m W b2 : invariant mass of the W + b or W −b system; • θ f 1 or θ f 2 : angles between fermion direction and opposite of the b orb quark in the W + or W − frame; • Φ f 1 or Φ f 2 : angle between the W + or W − decay plane and thetW + b or tW −b plane in the t ort quark frame; The decay of the H boson with angles θ * V and Φ * V is shown only for illustration, their distribution is flat for a spin-zero H boson production due to the lack of spin correlations between the production and decay processes. Their complete description is discussed in Ref. [61] in terms of the two equivalent angles θ * and Φ 1 . The grouping of the W + W − and bb systems, as opposed to W + b and W −b , is motivated by enhanced spin-correlation effects visible with the corresponding observables. The complete multidimensional distribution retains full information with either approach. Figure 5 shows the non-trivial angular distributions in the process pp → ttH corresponding to four scenarios of anomalous ttH couplings: pure scalar, pseudoscalar, and two mixed scenarios with f CP = 0.28 (corresponding to the equal cross-section of scalar and pseudoscalar processes) and different phases. Most angular observables exhibit a clear difference between the scalar and pseudoscalar processes. Only observables appearing in sequential decay of the top quarks are sensitive to φ CP . As noted earlier, these observables together with the boost of the ttH system are equivalent to other observables defined in the laboratory frame, as shown in Appendix A, but provide complete kinematics required as input to the matrix element tools and emphasize particular features in the process.
The description of observables Ω in the other processes pp → tqH and bbH follows by analogy, with only a subset of observables available due to lack of sequential decay of at least one associated quark.

B. Matrix element likelihood approach
With the kinematics of a process reflected in the complete set of observables Ω, one could in principle analyze the data in this multi-dimensional space. However, this often becomes impractical with a large number of observables, as illustrated in Section IV A, when parameterization of probabilities and detector effects in such a multi-dimensional space becomes difficult. Reducing the number of observables is a possible approach, but essential information may be lost. Machine learning techniques can approximate optimal functions that can depend on a large number of inputs, but those also require training and therefore perform only as well as training goes. These techniques are typically targeted to discriminate between certain categories of events and are not optimal for dealing with quantum mechanical interference effects which become essential in physics processes, though for possible solutions see Ref. [92].
The matrix element techniques are the methods based on multivariate per-event likelihoods prepared using a phenomenological calculation for the process of interest. They may employ the same calculation as used in the Monte Carlo event generators or may be reformulated to represent analytical distributions of observables of interest, such as Ω. Such matrix elements, if used properly, are guaranteed to retain full information about the event. The difficulty in using matrix element methods often comes from non-trivial detector effects which alter multivariate likelihoods. This issue is greatly simplified by utilizing ratios of the matrix elements in which certain detector effects cancel, most importantly variation of reconstruction efficiency as a function of observables. Resolution effects may be introduced with transfer functions, or neglected when their effect on performance is small. Missing degrees of freedom, such as neutrinos, may be either constrained from the global event information, as we illustrate below, or integrated out in the matrix element calculation. In the end, either machine learning or matrix element techniques could be used in the analysis of the data, but in either approach, it is ultimately the matrix elements which guide us in maximizing the amount of information, as they are also used in machine training through Monte Carlo.
The basic idea of the MELA technique is to project kinematics on the minimal set of discriminants calculated as ratios of the matrix elements. It has already proved to be a powerful tool for the H boson discovery and characterization during Run I of the LHC as applied primarily to the H boson coupling to the vector bosons. For a simple discrimination of two hypotheses, the Neyman-Pearson lemma [93] guarantees that the ratio of probabilities P for the two hypotheses provides an optimal discrimination power. However, for a continuous set of hypotheses with an arbitrary quantum-mechanical mixture several discriminants are required for an optimal measurement of their relative contributions. For example, probability for interference of two contributions could be presented as where P int and P ⊥ int describe quantum mechanical interference of J P = 0 + and 0 − terms. One could apply the Neyman-Pearson lemma to each pair of points in the parameter space of (f CP , φ CP ), but this would require continuous, and therefore infinite, set of probability ratios. However, an equivalent information is contained in a linear combination of only three probability ratios, which can be treated as three independent observables. For H boson physics at proton or lepton colliders, such discriminants are introduced in Ref. [63] as which become the optimal discriminants for the process with four contributions in Eq. (6). In Eq. (6), x i is a set of observables describing the process, which may be Ω when calculating the discriminants, or may be discriminants themselves when performing the analysis later. The number of discriminants can also be reduced by dropping D ⊥ CP assuming sin φ CP = 0, which is the case for real κ f andκ f in Eq. (2). On the other hand, with additional contributing amplitudes the number of observables grows. For example in the presence of background the D bkg discriminant is introduced which can also be supplemented by the interference discriminant if there is quantum mechanical interference between the signal and background processes. The corresponding two discriminants are defined as Calculating a discriminant analogous to D bkg for the pseudoscalar signal hypothesis is not necessary as a combination of Eqs. (7) and (10) carries the needed information. The number of discriminants grows with the number of free components in the model; for example the background may interfere with different signal components and those may require different observables. However, typically there is a limited set of interference discriminants which become of practical interest, as we illustrate below. The probabilities P in Eqs. (7)(8)(9)(10)(11) are the physical cross sections given by the product of parton distribution functions convoluted with the partonic cross sections that are proportional to the squared matrix elements. The latter depend on the full event kinematics as measured in the experiment or simulated by a Monte Carlo generator. They are computed at LO and do not include detector effects. However, as we illustrate in the following studies, they remain nearly optimal even after higher order or detector effects are introduced. The probabilities P in Eq. (6) may be treated as templates of the limited number of optimal discriminants when the analysis is performed. These templates are obtained from numerical simulation of the processes accounting for parton showering and detector effects. In the following analysis we limit the maximum number of discriminants to three, which we find to be both practical and close to optimal.
The complete set of optimal discriminants in Eqs. (7)(8)(9)(10)(11) was introduced earlier in experimental analysis of HV V processes with LHC data by the CMS [2,3,52,53,55,57,91] and ATLAS [5,59,60] experiments, and phenomenological studies supporting this development [61][62][63]. For example, it was shown that the complete set x i = {D 0− , D CP , D bkg } was optimal for the measurement of f a3 , a parameter equivalent to f CP , for the real HV V couplings. A subset of equivalent observables was also introduced independently in earlier work on different topics [94][95][96]. Here we apply this formalism to the measurement of the H boson anomalous couplings to the heavy flavor fermions for the first time.

C. Application to the ttH process
The large number of observables Ω defined in Section IV A for the ttH process can be compressed in a compact form with only three discriminants x i = {D 0− , D CP , D ⊥ CP } as defined in Eqs. (7)(8)(9), which is sensitive to the measurement of anomalous ttH couplings. The distributions for the discriminants are shown in Fig. 6 for J P = 0 + , 0 − , and mixed states. The nominal studies presented here are based on LO in QCD calculations. Variations due to NLO matrix elements, PDFs, and QCD scale uncertainties are also shown in Fig. 6 and are discussed in more detail in Section V.
As one can see from both the discriminant definitions and distributions in Fig. 6, the D 0− is sensitive to the relative size of CP contributions, while D CP and D ⊥ CP are sensitive to CP mixture leading to forward-backward asymmetry in the presence of both CP amplitudes for the real and complex ratio of couplings κ f /κ f , respectively. It is interesting to observe that the asymmetry is not strongly pronounced in the case of real couplings even when using the full top decay chain information. The asymmetry is more pronounced in the case of complex couplings, as seen in the D ⊥ CP distribution, which can be traced to the Φ f 1 distribution in Fig. 5. The asymmetry in both D CP and D ⊥ CP disappears when top decay information is not used in the matrix element, which reflects the fact that spin correlations in the tt system decay are essential for observing effects sensitive to CP mixture.
At the moment of discovery of the ttH process, precision will be limited by statistics and the D 0− discriminant will provide the most information about the CP components in the process. As smaller anomalous contributions get tested, the importance of the interference discriminant will grow, but ultimately the full information is contained in the complete set of discriminants.

D. Application to the bbH process
The bbH and ttH processes are very similar with the main difference being the heavy quark mass which, in fact, has a significant impact on the sensitivity of kinematic shapes to the Hff couplings. This is because shape sensitivity arises from the mixing of left-and right-handed helicities at the matrix element level. Therefore, this effect becomes proportional to the mass of the associated quark and becomes essentially invisible in the bbH process. In Fig. 7, we plot the angular distributions as well as the matrix element discriminant for the bbH process, analogous to ttH process distributions shown in Figs. 5 and 6. The different CP states have almost identical distributions, as follows from the helicity flip suppression discussed above. Therefore, we conclude that it will be very challenging to probe the CP nature of Hbb coupling through shape analyses in the bbH production mode.

E. Application to the tqH process
The tqH production process features both fermion and vector boson couplings of the H boson, as shown in Fig. 2. Interference between the Hff and HV V -induced diagrams in Fig. 2 is destructive in the SM, but any deviation in either size or sign of either contribution could lead to a significant change in observations. Therefore, in this paper we illustrate the approach where two parameters of interest are determined: relative size of the Hff and HV V contributions, including their relative phase, and the relative size of the anomalous Hff coupling. In this context, contributions from the HV V process could be considered as background and we consider only SM-like HV V coupling with a 1 = 0 in Eq. (4), while the signal process with the Hff coupling is allowed to have an arbitrary anomalous contribution. Therefore, the three discriminants D 0− , D bkg , and D int bkg as defined in Eqs. (7,10,11) provide the most relevant information for this analysis. Their distributions are shown in Fig. 8.
There is a clear difference between distributions for the alternative hypotheses, such as between J P = 0 + and 0− Hff -induced signal in D 0− , or between J P = 0 + signal and HV V -induced process in D bkg . It is important to stress that destructive interference between the J P = 0 + signal and HV V -induced processes with SM couplings leads to distributions which are very different from the direct sum of the two distributions, as shown in Fig. 8. In particular, the D int bkg discriminant shape is significantly distorted due to effect of interference, while the other two discriminants also exhibit sizable differences as well. This feature leads to strong separation power between different hypotheses even with small number of events in analysis, as we show below.

F. Application to the H → τ + τ − process
In the H → τ + τ − process, it is possible to define the full sequential decay kinematics and construct the matrix elements using information about all final state particles. This is illustrated for the process H → τ + τ − → ( νν)( νν) with the D 0− discriminant in Fig. 9. Even though there is a strong separation power between the J P = 0 + and 0 − models in this case, there is little practical application because reconstruction of the four neutrinos is not possible. Therefore, only limited information can be retained in reconstructed observables and we use the matrix elements for the purpose of MC re-weighting techniques below. In the case of hadronic τ decay, we provide the matrix element for the H → τ + τ − → (Xν)(X ν) process, where X could be any hadronic particle decayed from τ , e.g. π, ρ, a 1 . Figure 10 shows the D 0− discriminant constructed using this matrix element, in a hadronic final state. The events are decayed through TAUOLA, including hadronic formfactors for particle hadronization, for the J P = 0 + and 0 − models. In addition, the J P = 0 + events are re-weighted to the 0 − model using the MELA weights, which allow us to create any model with arbitrary anomalous couplings. The D 0− discriminant can be compared to other observables proposed for analysis of the H → τ + τ − decay, for example Φ CP [36], defined as q τ − and q X,X are the 3-momentum of τ − and X or X in the H boson rest frame. The two observables, D 0− and Φ CP , carry similar information for analysis of anomalous couplings, but the D 0− discriminant is somewhat more powerful.

V. NLO QCD STUDY OF KINEMATIC DISCRIMINANTS
Let us now discuss the effects of higher order QCD corrections on the modeling and performance of anomalous coupling discrimination. As described in Section III, we compute the NLO QCD predictions for pp → ttH production followed by the spin-correlated top quark decays at NLO QCD in the narrow-width approximation. Neglecting QCD corrections in the description of the pp → tt + H process constitutes the dominant theoretical uncertainty on its cross section. We find that residual scale uncertainty on the total cross section is reduced from 21 % to 9 % when going from LO to NLO in QCD. The corrections on shapes of basic kinematic distributions are up to ±10 %.
In the earlier work [63], we studied the impact of NLO in QCD effects in production on the anomalous coupling discrimination in decay H → V V . However, the production and decay processes carry no spin correlation and additional radiation from the production stage is largely decoupled from the color-neutral H → V V decay dynamics. Hence, it is straightforward to use LO matrix elements to characterize HV V couplings, even in the presence of initialstate radiation. This is in contrast to the pp → tt (→ bbW W ) + H(→ V V /ff ) process where initial and final state particles radiate color charges and the top quarks exhibit spin correlations, all of which affect the studied dynamics.
A fully consistent extension of the matrix element method beyond the LO requires both event generation and matrix element discriminants at higher orders. The main complication stems from kinematic configurations where hadronic activity is clustered by a jet algorithm. Commonly-used jet algorithms combine soft and collinear radiation in subsequent 2 → 1 recombination steps. Hence, the resulting jet either acquires some invariant mass which does not correspond to the fundamental parton mass, or the jet violates global momentum conservation. This feature prohibits the use of jet momenta in a LO matrix element, which has on-shellness (of quarks and gluons) and momentum conservation built-in from first principles. A systematic solution to this issue at NLO QCD is part of active research and first elegant solutions have been presented in Refs. [97][98][99][100], where modified jet algorithms are proposed to map resolved and unresolved parton configurations onto their proper matrix elements. These approaches have promising prospects for future measurements, but they require the use of new jet algorithms that are currently not used in experimental analyses. Moreover, only solutions for either colorless final states or colorless initial states have been presented in the literature. A fully developed application to e.g. top quark pair production at NLO QCD is not yet available. A variation away from the exact NLO treatment has also been presented in Ref. [101], where additional radiation is included through a parton shower approximation. This approach allows to include multiple emissions and has been applied to Higgs boson physics in Refs. [102,103].
In this paper, we take a pragmatic and more simplistic approach. We retain leading order matrix elements in the discriminants of Eqs. (7)(8)(9)(10)(11) and probe them with events from leading and next-to-leading order simulation, and also compare those to variations due to PDFs, QCD scales, and parton showering. The mismatch between the LO discriminants and NLO simulation does not formally allow us to claim optimal discrimination power by virtue of the Neyman-Pearson lemma, where constructed likelihoods should be interpreted as fundamental probabilities. However, we demonstrate that NLO corrections to the shapes of kinematic distributions in the pp → ttH process are small and sometimes indistinguishable when compared to other associated uncertainties. Therefore, the LO discriminants D maintain their discriminating power beyond the well-defined leading-order, and we can continue to use them as robust and powerful tools for anomalous coupling studies.
In Fig. 6, we compare the impact of LO versus NLO events probing the LO discriminants D. The solid histograms show the distribution for LO events, whereas the hashed bands indicate the shift due to NLO corrections. We note that the general shapes of the various distributions are maintained and only minimally distorted. The separation power between the extreme J P = 0 + and 0 − hypotheses is largely unaffected by the presence of higher order corrections. The most powerful discriminating observable D 0− receives very small corrections in range within the bulk of the distributions, as shown in the lower pane of Fig. 6. Moreover, most of this correction appears already with the PDF variations before NLO corrections at the matrix element level. Hence, the bulk of corrections that we observe stems only from different input parameters and PDFs. The width of the bands in all lower panes of Fig. 6 corresponds to scale variations by a factor of two around the central scale µ = m t + m H /2. Studies presented in Fig. 6 do not include parton showering. However, as we show in Section VI A and Fig. 11, inclusion of parton showering in LO simulation brings this simulation even closer to NLO modeling with parton showering.
We therefore conclude that discrimination power of the MELA approach is guaranteed even when higher-order corrections are considered in the pp → ttH process and additional jets are present in the event sample. Soft and collinear radiation, which generates massive jet momenta, can be handled in the matrix element approach and does not spoil the discrimination power. These higher-order effects are within the uncertainties of the PDF, scale choice, and parton showering.

VI. APPLICATION TO CP PARITY MEASUREMENTS IN ttH, tqH, AND bbH
In this section, we estimate the potential for CP measurements in the ttH, tqH, and bbH processes on LHC with 300 fb −1 of proton-proton collision data collected at 13 TeV center-of-mass energy. This is the integrated luminosity expected by the end of Run III of LHC in about seven years. Projections to other luminosity scenarios are usually trivial extensions as long as uncertainties remain limited by statistics. While there is a strong evidence for the ttH production in Run I of LHC, none of these processes have been firmly established yet. However, we rely on experimental studies of these processes in Refs. [39][40][41][42][43] for realistic event reconstruction projections.
As the first observation, following Section IV D, we conclude that it will not be possible to measure CP in the bbH production process in the LHC program. For the ttH and tqH processes, we consider the H → γγ decay mode to tag the H boson, as a clean signature with sizable branching fraction. We also consider the H → ZZ → 4 final state in the ttH process for comparison, but its contribution is small due to the small branching fraction. We use the hadronic decay of one top quark final state so that the full kinematics can be reconstructed. In the ttH case, the other top quark is reconstructed in the leptonic channel. Inclusion of other final states of either H boson or top would only enhance expected precision, but the decays we consider are representative of the typical analyses of these processes.
In this study, the ttH, bbH, and tqH processes with SM or anomalous couplings are generated with the JHU generator. The only non-negligible background that we need to consider is SM ttγγ production as a background to the ttH study with the photon decay of the H boson, which is simulated with MadGraph. The MC samples are interfaced to PYTHIA8 for parton shower and hadronization. In order to model detector effects, lepton and photon p T are smeared with 1% and 4% resolution. The jets are reconstructed in a cone of R = 0.5 using anti-k T algorithm and their energy is smeared by 20%.
The event selection criteria follow those of the LHC analyses [40]. We require the leptons, photons, and jets to have p T > 5, 10, and 30 GeV, and |η| < 2.4, 2.4, and 4.7, respectively. Jets within ∆R < 0.2 of the leptons or photons are removed. In the ttH analysis, an event should have at least four jets and a b-tagged jet. The b-tagging efficiency (62%) and fake rate for the light quark jets (6%) follow experimental study [40]. To fully reconstruct the semileptonic decay of the tt system, we use the constraint fit from Ref. [104]. The four-momenta of four jets, MET, and one lepton are used in the kinematic fit with the masses of the top quarks and the W bosons as constraints. If more than four jets are reconstructed, the combination that gives the best χ 2 is selected. The four-momenta of all decay products of the tt system are obtained from this fit and are used in the further analysis. In the tqH analysis, exactly four jets and a b-tagged jet are required in order to remove hadronic ttH events. The combination of three jets with the mass closest to the top is treated as the top decay product in this process. The required number of reconstructed leptons and photons depends on the studied final state. If required, the leading photon should have p T > 33 GeV and p T /m γγ > 0.5. In the H → 4 channel, two pairs of opposite sign and same flavor leptons should have invariant mass greater than 40 and 12 GeV. The invariant mass of the H boson candidate is required to be between 100 and 140 GeV.
In the case of the tqH process with H → γγ, the main other contribution is cross-feed from the ttH process with the same H boson decay. The ttH process with the 4 decay of the H boson has negligible background, while with the γγ decay the dominant background is the SM ttγγ production. The expected number of events of signal and background events at 300 fb −1 is shown in Table II. We would like to note that these expected yields are quoted for the SM scenario where destructive interference between the Hff and HV V -induced tqH processes leads to a small number of expected events. However, this interference may become constructive with the non-SM couplings. The cross section for ttγγ processes suggests background yield to be smaller than the signal. However, the LHC studies with data-driven methods suggest larger background [40]. Therefore, we conservatively set the ttγγ background yield to be twice the signal in the invariant mass window specified above. The analysis of the ttH process uses the D 0− discriminant, where decay of the top quarks is not considered in the matrix element. Consideration of the top quark decays is important in the calculation of the D CP or D ⊥ CP discriminants, but only when the up and down flavors of the quarks in the decay chain is known. The latter is difficult to determine with the jet reconstruction techniques and therefore the CP discriminants are not used in this analysis. In the H → γγ channel, we use the invariant mass m γγ to separate the signal and background. Figure 11 shows the D 0− distribution in the H → γγ channel, where the J P = 0 + , 0 − and background distributions are shown.
In Fig. 11, simulation of the 0 + process is also shown with the POWHEG generator at NLO in QCD. In all cases, parton showering is performed with PYTHIA8. Similar to the study presented in Section V, the NLO QCD effects are found to have a small effect on accuracy of D 0− simulation, especially after parton showering is included in simulation. Any residual effects are consistent with systematics also arising from PDF and QCD scale variations.
The expected precision of the f CP measurement in the ttH process with both H → γγ and H → 4 decays, and their combination, is shown in Fig. 11 for integrated luminosity of 300 fb −1 . The maximum likelihood fit is based on the probability density functions following Eq. (6) parameterized with template distributions filled with generated events as discussed above. About 3σ exclusion of the pure pseudoscalar state is expected in such a scenario, which is comparable to the current precision with the HV V measurements, but provides a fundamentally different approach through fermion couplings. Scenarios with a sizable CP mixture, |f CP cos φ CP | > ∼ 0.8, are excluded at 2σ.

B. Study of the tqH process
The analysis of the tqH process uses the D 0− , D bkg , D int bkg discriminants, shown in Fig. 12. In this study, the Hff -induced process is considered as signal and the HV V -induced process is considered as background. Similar to the ttH study, the decay of the top quarks is not considered in the matrix element and the D CP discriminants provide little information and therefore are not used. There is a sizable contribution of the ttH events misreconstructed as tqH and they carry information on f CP . The above observables provide sufficient information to differentiate between CP components of the ttH process as well. All event contributions in this study can be parameterized with three couplings, κ,κ, and a 1 , which are assumed to be real, as where L is integrated luminosity and σ is the product of cross section and reconstruction efficiency of a particular process corresponding to the unit value of the coupling (κ,κ, or a 1 ). The interference cross section can be negative, as it is for interference between the κ and a 1 terms in the tqH process. In the tqH process, we express κ,κ, and a 1 in terms of two effective cross section fractions with phases and the overall normalization as follows fκ =κ 2 σ tqH 0− (a 2 1 σ tqH bkg + κ 2 σ tqH 0+ +κ 2 σ tqH 0− ) , φκ = arg(κ/a 1 ) = 0 or π .
In the SM, f κ = 0.46, fκ = 0, and φ κ = 0. The ratios of cross section is σ tqH 0+ /σ tqH bkg = 0.86. The maximum likelihood fit, similar to the ttH analysis, uses a 3D template approach of three observables D 0− , D bkg , D int bkg , and with f κ , fκ, and total event yield as free parameters. The expected precision of the fit is shown in Fig. 13 (left plot). This approach allows simultaneous measurement of both the relative fraction of HV V and Hff induced processes and of the anomalous contribution in the Hff coupling, with proper accounting for all interference effects. This measurement can be reduced either to the measurement of f κ with the constraint fκ = f CP = 0 (middle plot), or to the measurement of f CP (right plot). Precision on the Hff couplings is driven by both tqH and ttH processes in this analysis, as illustrated with the likelihood scans separated for the two samples of events in the right plot of Fig. 13. More than 3σ exclusion of the pure pseudoscalar state is expected in such a scenario, which is a measurement independent from that discussed in Section VI A since ttH events have little overlap. It is important to note that in this scenario it will be possible to determine the relative sign of the Hff and HV V -induced contributions and exclude both extreme scenarios of either pure Hff or pure HV V processes, assuming events follow SM expectation. , and D int bkg (right) discriminants in the tqH study where the Hff -induced process is considered as signal and the HV V -induced process is considered as background. The following three contributions are considered: J P = 0 + signal (red crosses), 0 − signal (blue triangles), HV V -induced process as background (black circles). Also shown are distributions of mis-reconstructed ttH signal with J P = 0 + (magenta diamonds) and 0 − signal (green squares). All distributions appear after simulation and reconstruction discussed in text.

VII. SUMMARY AND CONCLUSIONS
We have developed the Monte Carlo simulation and matrix element analysis tools and investigated prospects for measurement of anomalous interactions in the H boson production in association with top or bottom quarks at the LHC, as well as its decay in two tau leptons. The study is based on the JHU generator framework and the matrix element MELA analysis technique. We find that it is difficult to measure anomalous couplings in the bbH process, while in both ttH and tqH analyses it is possible to have more than 3σ separation of the pseudoscalar hypothesis from the scalar with 300 fb −1 of integrated luminosity at LHC at 13 TeV. It is also possible to separate the Hff and HV V processes and determine their relative phase in the tqH production, where in the SM the two processes interfere strongly and destructively. This feasibility study considers only representative decay channels of the top quark (hadronic decay of one top) and H boson (diboson decay) and inclusion of other final states would only enhance expected precision. Systematic uncertainties from QCD effects, such as PDF, scale, parton showering, and higher order corrections, are shown to be relatively small compared to expected statistical precision. The tools and techniques presented should facilitate measurements of SM and anomalous Hff couplings. ment and thank Roberto Covarelli, Chris Martin, Candice You for help with the generator validation. We are grateful to our co-authors of the JHU generator project for valuable contributions, and in particular to Ulascan Sarica and Heshy Roskes for the generator package support and Fabrizio Caola for providing useful comments on the manuscript. This research is partially supported by U.S. NSF under grant PHY-1404302 and by the German Federal Ministry for Education and Research (BMBF) under grant 05H15VKCCA. Calculations reported in this paper were performed on the Homewood High Performance Cluster (HHPC) of the Johns Hopkins University and the Maryland Advanced Research Computing Center (MARCC).
Appendix A: Supplemental information on kinematics Figure 14 shows the kinematic observables defined in the laboratory frame in the SM process gg and qq → ttH, corresponding to four scenarios of anomalous ttH couplings. These observables can be derived from those shown in Fig. 5 and defined in Section IV A. The kinematic distributions in the process gg and qq → ttH defined in the laboratory frame: top quark rapidity (Y (t)), transverse momentum of the top quark (pT (t)), H boson rapidity (Y (t)), transverse momentum of the H boson (pH (t)), ttH system rapidity (Y (ttH)), pseudorapidity difference between the two down type fermions decayed from top and anti-top (∆η ll ) and between two bottom quarks (∆η bb ), cos θ ll between the two down type fermions, and cos θ bb between the two bottom quarks. Four scenarios of anomalous ttH couplings are shown: fCP = 0 (SM 0 + , black circles), fCP = 1 (pseudoscalar 0 − , red crosses), fCP = 0.28 with φCP = 0 (blue triangles) and φCP = π/2 (green diamonds). The LHC pp energy of 13 TeV and H boson mass of 125 GeV are used in simulation.