Probing the $CP$ nature of the top quark Yukawa at hadron colliders

We analyze the prospects of probing the $CP$-odd $i \tilde \kappa \bar t \gamma^5 t h$ interaction at the LHC and its projected upgrades, the high-luminosity and high-energy LHC, directly using associated on-shell Higgs boson and top quark or top quark pair production. To this end we first construct a $CP$-odd observable based on top quark polarization in $W b \to t h$ scattering with optimal linear sensitivity to $\tilde \kappa$. For the corresponding hadronic process $pp \to t h j$ we present a method of extracting the phase-space dependent weight function that allows to retain close to optimal sensitivity to $\tilde{\kappa}$. We project future sensitivity to the signal in $pp \to t(\to \ell \nu b )h(\to b\bar b) j$. Based on these insights we propose novel $CP$-odd observables for top quark pair production in association with the Higgs, $pp \to t \bar t h$, with semileptonically decaying tops and $h\to \bar b b$, that rely solely on measuring the momenta of leptons and $b$-jets from the decaying tops without having to distinguish the charge of the $b$-jets. Among the many possibilities we single out an observable that can potentially probe $\tilde \kappa \sim 0.1$ at the high-energy LHC with $2\sigma$ confidence.


I. INTRODUCTION
The coupling of the 125 GeV Higgs boson (h) to the top quark, which is the largest of the Standard Model (SM) couplings, is an important target for the LHC experiments. CP -violating h couplings are particularly interesting as any sign of CP violation in Higgs processes would constitute an indisputable New Physics (NP) signal. Existing data on Higgs production and decays is already precise enough to constrain any isolated modification of the top Yukawa to O(1) [1,2]. However, all existing measurements are based on CP -even observables with very limited sensitivity to CP -odd modifications of the top quark Yukawa. In principle, indirect collider bounds from Higgs decay and production (gg → h, h → γγ), and especially the low-energy bounds on electric dipole moments (EDMs) atoms and nuclei that target specifically CP -odd effects [1,3,4], are currently more constraining than direct collider probes. However, these constraints are subject to assumptions about other Higgs interactions, and in particular in the case of EDMs also other contributions unrelated to the Higgs.
A few existing proposals for LHC measurements of top quark pair production in association with the Higgs boson have studied manifestly CP -odd observables with on-shell t,t and h [1,[4][5][6][7][8][9][10]. In particular, the CP nature of the top-Higgs coupling in this case is reflected in the correlation between the spins of the tops, which can be reconstructed using the angular distributions of the top quark decay products. It turns out however that the resulting effects are typically almost prohibitively difficult to measure at the LHC due to limitations of simultaneous top quark reconstruction, as well as their spin and charge identification. An alternative is offered by the single top production with associated Higgs based on the hard process bW → th and observed as pp → tHj. Owing to simpler kinematics in this case the measurement of top quark polarization is accessible through the angular distribution of lepton in the top rest frame. Several existing studies in this direction have already proposed top quark polarization related observables [1,[11][12][13][14]. Yet while the literature abounds with proposals of CP -sensitive measurements both in tth and th channels, there has been no study to systematically search for and construct observables with optimal sensitivity to CP-odd top Yukawa under realistic conditions at hadron colliders.
In this paper we address this challenge by identifying observables with optimal sensitivity to a single CP -odd parameter in both th and tth associated production at the LHC, which can be realistically measured and exhibit close to optimal sensitivity to CP -odd interactions between the Higgs boson and the top quark. The proposed observables in th are based on optimization of top-spin correlations previously studied in tt production [15]. In the case oftth this procedure becomes intractable in practice and our construction relies instead on CP -and P -symmetry arguments. To set the stage we write the effective top quark -Higgs boson interaction as where y t = √ 2m t /v is the top quark Yukawa in the SM, while real dimensionless quantities κ,κ parametrize departures from the SM (at κ = 1,κ = 0). In the context of the SM Lagrangian complemented by dimension-6 effective interactions,κ is generated from the term |H| 2QH u R which is misaligned relative to the SM Yukawa [16]. Clearly, any indication of non-vanishingκ would be an indisputable sign of NP. Our goal is to construct optimized and practically measurable observables which probe the CP -odd parameterκ directly.
The rest of the paper is structured as follows. In Sec. II we study optimized top spin observables in associated single top and Higgs boson production, both in idealized partonic W b → th scattering, which is tractable analitically, as well as in more realistic simulations of thj production and reconstruction in proton collisions at pp colliders. Sec. III contains the analysis of CP violating observables built from accessible momenta in tth production, both at Monte Carlo (MC) level and after including the background, reconstruction, and detector effects. Finally, our main conclusions are summarized in Sec. IV.
We begin by studying the effects ofκ on top spin observables in the idealized case of W (p)b(p ) → t(k)h(k ) scattering, where the complete polarized scattering amplitudes can be found in a compact analytic form. This process can actually be connected to a more realistic pp → thj production in the high energy limit, where the W and b quark mass effects are negligible and the collinear emission of both initial state 'partons' can be described by the corresponding parton distribution functions. Three diagrams contribute to such parton level Higgs-top production in the SM, shown in Fig. 1. Neglecting furthermore the mass (and thus the corresponding Yukawa coupling) of the bottom quark, we consider only the first two of the diagrams in Fig. 1. The formalism presented here is based on Refs. [15,17,18]. First, we introduce the spin projection operator where s µ is a top spin four-vector, defined in a general frame as Vector k is the top quark momentum andŝ is an arbitrary unit vector. The physical significance ofŝ is revealed if we make a rotation-free boost 1 to the top rest frame where we find s * µ = (0,ŝ). Therefore s 2 = −1, s · k = 0, andŝ corresponds to the polarization of the top quark in its rest frame. Projection onto a well defined polarization of the top quark is achieved by inserting the operator (2) at the amplitude level and leads to the following relation at the cross-section level: Thus the cross-section is linear in s µ where b µ contains all information about the polarization of the top in the process. The parton level cross section can be writen as where Φ i is the initial state flux normalization and dΓ th is the th phase space volume. On the other hand, in the top rest frame it is convenient to introduce the spin density matrix as such that the unpolarized cross section is proportional to |M| 2 = Tr[ρ] = 2A. Here σ are the Pauli matrices. In the density matrix formalism, the expectation value of a generic operator is obtained as In particular, the polarized cross section alongŝ is obtained as the expectation value of the projector: One can determine the rest-frame coefficients A, B i from a, b µ by comparing the expressions for polarized |M| 2 , expressed via Eq. (5) and Eq. (9). The result of this matching are explicit expressions: The rest-frame polarization of the top quark along a vectorŝ is given by the expectation value of Oŝ =ŝ· σ 2 , which can be probed by measuring the angular distribution of the lepton in the top decay 2 thus allowing for experimental extraction of B i coefficients [19]: Here θ is an angle between the lepton and the polarization axisŝ in the top rest frame. The above construction shows that the vectorŝ is an arbitrary unit vector defined in the laboratory frame. A particular choiceŝ = k/|k| implies that Oŝ measures the top quark helicity. Another natural choice forŝ is the W momentump, also known as the beam basis, which has to be redefined in pp collisions where the W momentum has a discrete ambiguity. Experimentally one has to reconstruct the top quark rest frame in order to be able to trace the angular distribution of the lepton with respect to the chosen s and gain access to the coefficients B i . In the following we will optimize the choice ofŝ such that the sensitivity to the CP -violating parameterκ is maximized.
In the W b center-of-mass frame we can define the W and t momenta aŝ where θ is the angle between the direction of the top quark and the W boson. We have set the azimuthal angle φ = 0 without loss of generality. The polarization vector components B i in this case depend on x = cos θ, and we have found that in the coordinate system in Eq. (13) the analytical expression for B 2 (x) is linear inκ, i.e., B 2 (x) = β 2 (x)κ, whereas B 1,3 do not contain linearκ termsEffectively this means that we should choose the vectorŝ to be be orthogonal to the plane spanned by the W and t momenta in order to probeκ with linear sensitivity. Similar results have been found in Ref. [13]. We fixŝ =p ×k/|p ×k|.
In this case the interesting experimental quantity is the following two-fold differential cross-section where we have approximated the intermediate top quark as a narrow resonance and Σ(x,ŝ) = dσ/dx(W b → t (ŝ) h) is the differential production cross section for the top quarks polarized in theŝ direction. Using Eq. (9) and insertingŝ we have Σ(x, ±ŝ) = Φ i (A(x) ±κβ 2 (x)). Thus we can write Eq. (14) as Treatingκ as a small perturbation we can integrate the distribution in Eq. (15) with a phase-space dependent function f that would maximize statistical sensitivity of the integral toκ. It has been shown in Refs. [20,21] that such an optimal function should be the ratio of theκ-perturbation to the unperturbed distribution, in our case f (x, cos θ ) = β2(x) A(x) cos θ . The optimal observable is thus where θ is the angle betweenŝ and the lepton momentum in the top center-of mass-frame, as defined in the preceding paragraph. The prediction scales as β 2 2 , where we have integrated over cos θ and left the bounds for x = cos θ unspecified. The function β 2 (x) is plotted in Fig. 2.
To carry over the presented formalism to the realistic case of pp collisions, we have to adapt the beam axis by referring only to experimentally accessible momenta. Using the reconstructed top momentum k as a reference, we define the positive z-direction as the parallel top quark momentum projectionk . The top quark is then always in the positive hemisphere,x = cosθ ≥ 0, whereθ is the angle between k and k . The polarization direction with linearκ sensitivity now becomesŝ =k ×k ⊥ upon which we now measure the lepton angleθ . The cross-section distributions inx and x are related via whereÃ The cos θ is flipped in the second term since forx = −x the polarization vectorŝ =k ×k ⊥ flips the direction compared to the previous definition,ŝ ∼ p × k. The optimal observable in this case is finallỹ . However in general theÕ W b→th opt.
is expected to result in a weaker statistical significance due to our inability of determining the direction of the top quark with respect to the initial W . Fig. (2) shows that β 2 (x) matches reasonably well the ideal scenario, β 2 (x) = −β 2 (−x), for representative values of the center-of-mass energy √ s.

B. Hadronic process pp → thj
Here we demonstrate the procedure of measuring the optimal observable in the case of pp collisions, but still neglecting reconstruction efficiencies and backgrounds. The parton level observable defined in Eq. (20) can be adapted to this case with an additional integration over the parton distribution functions (PDFs). Since the hadronic cross section is a convolution of partonic cross sections it can be split into aκ-independent piece and the small perturbation proportional toκ, similar to the partonic cross section in Eq. (18). Assuming that the Higgs decays into visible states, the missing p T is only due to the neutrino originating from the top decay. Thus we can reconstruct the top quark momentum and kinematic quantities of Eq. (18). Thus, for hadronic collisions one can express the cross section as and weigh the events with the optimal f opt. ∼ cosθ B 2 /A. We use the MC event generator MadGraph5 [22,23] together with the Higgs Characterisation UFO model [24,25] (for analysis including NLO QCD effects see also Refs. [6,12], for NNLL EW effects [26]) to incorporate the κ andκ couplings in the simulation of the pp → t(→ b ν)hj signal. The procedure of extracting the weight function B 2 /A from MC simulations and using it to produce the optimal observable goes as follows: 1. Choose the bins forx betweenx min ≥ 0 andx max ≤ 1.
2. Fixκ and extract from the MC simulation the mean cosθ in each of thex bins. The obtained value corresponds to weight 1 3 B 2 /A in this bin, see Eq. (21).
3. Use this information to weigh experimental events bin-by-bin with f opt. ∼ cosθ B 2 /A. The normalization of f opt. is fixed by the requirement dxB 2 /A = 1.
This optimization procedure is independent of theκ value. The resulting optimal weight B 2 /A is shown in Fig. 3, where we compare it for different final states (thj orthj) and collision energies (14 or 27 TeV).
We have also extracted the weight function including NLO QCD simulations to estimate the systematic uncertainty and found that the difference is within 10% of the LO extraction. Contrary to the presented approach one can accessκ by measuring cos θ , which corresponds to the case where the weight is independent ofx, i.e. B 2 /A = 1. Fig. 4 shows improvement of the significance when the optimal weight function is applied on simulated signal events without showering or reconstruction effects at 14 TeV.  In order to make contact with experiments, we now present realistic bounds in the (κ,κ) plane. For this we include the effects of parton showering, the detector response and background proccesses. We use MadGraph5 to generate events at leading order (LO) in QCD for the signal process pp → t(→ b ν)h(→ bb)j plus the conjugate process witht at 14 TeV High-Luminosity LHC (HL-LHC) and 27 TeV High-Energy LHC (HE-LHC) center-of-mass energies 3 . Event generation is performed for multiple values of (κ,κ). The parton level events are subsequently showered and hadronizied with Pythia8 [27], and jets are clustered with the anti-k T algorithm using FastJet [28]. For detector simulation and final state object reconstructions (e.g. lepton isolation and b-tagging) we use Delphes [29] with the default ATLAS parameters in delphes card ATLAS.tcl from version Delphes-3.3.3. The dominant background process in this analysis is tt plus jets. We include this background by generating pp → tt samples decayed into the semi-leptonic channel, produced in association with 0, 1 and 2 hard jets. In order to correctly model the hard jets, we merge the matrix element computations with the MC shower using the MLM [30] jet matching algorithm. For the event selection we demand the following basic requirements: • Exactly 3 b-tagged jets with |η(b)| < 5 and p T (b) > 20 GeV, • One additional (non-tagged) light jet exclusively in the forward direction with 2 < |η(j)| < 5 and p T (j) > 20 GeV, • One isolated light lepton ± = e ± , µ ± with |η( )| < 2.5 and p T ( ) > 10 GeV.
In addition, we further select events with one reconstructed Higgs and one reconstructed top quark as follows: first, we calculate the three possible invariant masses from the three reconstructed b-jets (m bb ) and only keep the event if at least one bb pair satisfies |m H − m bb | < 15 GeV. For such events, we select as the Higgs decay candidate h → bb for the pair of b-jets with the invariant mass closest to the Higgs mass. The remaining non-Higgs b-jet is then assumed to come from the top-quark decay. Next, we reconstruct the top-quark by requiring that the combined invariant mass m blν of the remaining b-jet, the lepton, and the neutrino (also reconstructed by assuming it to be the unique source of missing energy in the event) to fall inside the mass window of the top-quark defined by m t ± 35 GeV. In order to further reject the tt backgrounds, events with a reconstructed Higgs and top are selected if the combined invariant mass of the b-jets originating from the Higgs and the light jet satisfies the cut m bbj > 280 GeV [31]. The final selection efficiency for the thj signal in the SM is 0.32% (0.23%), while for the background it is 0.008% (0.006%) at 14 TeV (27 TeV).
As we fully reconstruct the th system and have access to the lepton momentum from the top decay we have all the necessary information for measuring the optimized spin observable. We use the optimal weight function B 2 /A (Fig. 3) extracted from the MC simulations to construct a χ 2 with an appropriately weighted signal process. Our results for pp → thj generated in the SM are given by the 2σ exclusion limits (shaded blue) shown in Fig. 5 for the HE-LHC at a luminosity of 15ab −1 . We also present the limits (given by the dashed black contour) assuming a 2σ positive excess above the SM expectation corresponding to a measurement of the optimized spin observable of O opt. = 0.06 ± 0.03 whose size and error are statistics-driven. Because of the nature of our observable, the signed fluctuation gives rise to asymmetric limits in theκ direction. In the κ direction the bounds are also not symmetric as pp → thj production is sensitive to sgn(κ). In order to include effects of the backgrounds, the same statistical analysis was repeated but this time combining the signal and the tt background in the χ 2 fit. As a result we obtain a signal significance of S/ √ B ∼ 0.8 (3.2) at 14 TeV (27 TeV) with a luminosity of 3 ab −1 (15 ab −1 ). Notice that even with a large background rejection, the irreducible backgrounds are simply too large and the signal is completely diluted. We leave the possibility of further optimizing the cuts in order to reduce the backgrounds or including other Higgs decay channels for future works.

III. CP -ODD OBSERVABLES IN pp → tth
In this section we consider new CP obervables in the process pp → tth, with the two top quarks decaying leptonically. This process has been recently discovered by the LHC collaborations [32,33] (for the state of the art prediction of the differential distribution see e.g. Ref. [26]) The top quarks in this process are known to be unpolarized, independent of theκ value [1]. On the other hand, the correlation among top spins does contain information on the underlying κ andκ parameters. We do not follow this approach in order to avoid combinatorial difficulties when reconstructing both t andt rest frames. Instead, we focus on lab frame kinematic distributions in variables which are CP -and P -odd and are constructed from accessible final-state momenta [4].  Table I: Momenta with well-defined C and P eigenvalues. The b,b, + and − are the top decay products. The last column is a rank-2 tensor -a direct product of an axial and a polar vector.

A. Laboratory frame CP -odd observables
We denote the 3-momenta of the leptons and b-jets originating from t andt with p + , p − , p b and pb, respectively, and the Higgs 3-momentum with p h . The C and P transformation properties of six independent combinations of these momenta are given in Tab. I. We focus only on combinations that are nontrivial under C, P (i.e., we omit scalars products) and are accessible in a realistic experimental environment. For example we consider p b + pb, but not p b − pb as differentiating between b andb is difficult experimentally (see Refs. [34,35] for recent attempts in extracting the charge of the b-jet). The six combinations of momenta in Tab. I are taken as a basis for constructing P -and CP -odd variables ω. This is achieved by contracting (anti)symmetrically the momentum tensors such that the resulting ω is C even and P odd. i.e. a pseudoscalar. The resulting spectrum is then linear in the pseudoscalar ω with the coefficient in front linear inκ, analogous to expression (21). At leading order inκ we find: In Eq. (22) we have parameterized the phase space with the pseudoscalar variable ω, while all other variables are collectively denoted by x. Now we can again exctractκ with the statistically optimal weight function, which in this case is given by f opt. ∼ γ(x)/A(x), while the associated observable is Here N is the number of experimental events. In contrast to the extraction of f for the th process, here the extraction of γ(x)/A(x) turns out to be more complicated due to high-dimensionality of the phase space (4 variables for each pp → tth, t → b + ν,t →b −ν ). One could use a Monte Carlo event generator to obtain events following the distribution in Eq. (22) in order to extract γ(x)/A(x). However, since binning in all dimensions is not feasible it is better to formulate the task as a maximization problem to obtain the unknown weight function f (x; α). Given the N events (x (i) , ω (i) ), i = 1, . . . , N , generated with a non-zeroκ, the corresponding observable and its associated standard deviation are obtained as respectively. The goal is to find the set of parameters α of the function f (x; α), defined on phase space x and parameterized by α, that maximize the significance: The significance Sig(α) is independent of large enough N . The arguments of function f could be scalar products between final state momenta, whereas the functional form, controlled by parameters α, should be general enough. The obtained f , that was optimized using MC data can then be applied on a given experimental sample. In the following we will not pursue the optimal weight but will perform partial optimization along a single dimension of phase space.
In the following paragraph we introduce the CP -and P -odd variables. The simplest pseudoscalar is a mixed product of the form V · A, where V (A) denotes vector (axial vector), an example of which is (p b − pb) · (p − × p + ) presented already in Refs. [4,36] (see also observables proposed in Ref. [8]). In our case we do not wish to use p b − pb which leads us to an alternative mixed product that does not rely on separating b fromb experimentally: Once we allow for a more complicated pseudoscalar of the form (V · A) (V · V ) there are 13 possibilities, listed in Appendix A. Out of those and the mixed product in Eq. (26), one variable stands out as the most sensitive one: The pseudoscalar variable ω 6 is bounded 4 within the interval [−1, 1]. We have found that for the differential cross section d 2 σ/(dx dω 6 ) (see Eq. (22)) the ratio γ(x)/A(x), where x is an arbitrary kinematic variable, is approximately constant and does not oscillate in sign, which allows us to use a naïve weight function, f (x) = 1, without paying too much price for the cancellation between contributions from different regions of phase space. The observable we use is thus the average of ω 6 : The behaviour of O 6 as a function ofκ is shown in Fig. 6.
In addition to the analysis of O 6 presented here, we have also analyzed other observables related to the other pseudoscalar variables in App. A. For some of them we have used the optimization technique along a chosen dimension of phase space x, which drastically improved their sensitivity to κκ, however not a single one of them reached the sensitivity of O 6 . We note that all discussed observables can be improved in sensitivity by performing a full optimization using Eq. (25).

B. Limits in the (κ,κ) plane from pp → tth
We now demonstrate the capability of current and future colliders to measure the O 6 observable in tth production. With MadGraph5, we generate multiple event samples of pp → tth for different values of (κ,κ), followed by the decay chain (t → b + ν ,t →b −ν , h → bb) at 14 TeV and 27 TeV. The partonic events are then fed into Pythia8 for showering and hadronization and finally into Delphes for detector simulation with the default ATLAS card. We follow the same steps to generate the events of the main irreducible background pp → ttbb, (t → b + ν ,t →b −ν ). The basic event selection requirements for this analysis are: • 4 or more jets of any flavor with |η(j)| < 5 and p T (j) > 20 GeV.
• Of which, at least 3 are b-tagged.
• Exactly 2 oppositely charged light leptons with |η( )| < 2.5 and p T ( ) > 10 GeV. Furthermore, in order to identify the b-jets from the top-pair decays we count the number of tagged b-jets N b and perform the following selections: if N b ≥ 4, we compute the invariant masses m bb of all possible b-jet pairs and select the pair with invariant mass closest to the Higgs mass m h = 125 GeV. If the selected pair falls inside the Higgs mass window defined by m h ± 15 GeV we remove the pair from the list of b-jets and select from this list the highest p T b-jets as our candidate top quark decay b-jets. However if N b = 3 we compute all possible invariant masses m bj where j are non-b-tagged jets in the event. We select as the h → bb candidate the bj pair that minimizes |m h − m bj | and falls inside the Higgs mass window m h ± 15 GeV. The remaining two b-jets are taken as the candidate top quark decay b-jets.
The reconstruction efficiency of signal events using this approach is 5% (4%) and for background events it is 4.7% (3.8%) at 14 TeV (27 TeV). We construct a χ 2 for the combined signal and background events. In this case the signal-to-background ratio is much more favorable with S/ √ B ∼ 33 (118) at 14 TeV (27 TeV), so a joint analysis is possible. Results for the 2σ exclusion regions in the (κ,κ) plane are shown in Fig. 8 for different integrated luminosities at 14 TeV (left panel) and 27 TeV (right panel). These results show that the HL-LHC is somewhat limited by statistics while the HE-LHC gives a promising coverage of parameter space, in particular it is sensitive to CP -odd couplings of order O(0.1) at high luminosities. In Fig. 8 we provide the 2σ and 5σ exclusion regions of the HE-LHC at 15 ab −1 . In order to illustrate the sensitivity to the sign ofκ, we also provide the same exclusion limits in the left panel (right panel) of Fig. 9 in the scenario where the measured central value of the observable is O 6 = (1.2 ± 0.6) × 10 −3 (O 6 = (2.4 ± 1.2) × 10 −4 ) at 14 TeV (27 TeV), where the quoted fluctuations and standard deviations are estimated from the statistical error. These measurements, corresponding to a 2σ excess over the expected null value in the SM.

IV. SUMMARY AND CONCLUSIONS
In order to establish, directly and with minimal additional assumptions, the presence of a CP -odd component of the top quark Yukawa (κ), we have studied manifestly CP -odd observables in th and tth production at the LHC and its prospective upgrades.
For the thj final states we have relied on the possibility of reconstructing the t quark momentum and accessing the t polarization. We have identified a particular polarization direction which is perpendicular to the th plane, where the top polarization along this direction would undoubtedly point to the presence of the CP -odd couplingκ. We have presented a method for optimizing the phase space dependent weight and shown its sensitivity at the HL-and HE-LHC for the semileptonic top and h → bb mode. The handful of signal events offer discriminating power, sensitive to the sign ofκ, however the irreducible background due to tt+jets severely dilutes the sensitivity of the proposed observable.
On the other hand, tth production has a considerably larger cross section at LHC energies compared to thj, while suffering more moderately from irreducible backgrounds. Due to the complexity of the final state kinematics with multiple undetected particles we have in this case proposed variables that only depend on the lab-frame accessible momenta and are manifestly P -and CP -odd. We have identified a single triple product variable that does not rely on b-jet charge determination. Finally, among the possible pseudoscalar variables constructed as products of five lab-frame momenta, we have singled out the most sensitive one, O 6 of Eq. (27), the sensitivity of which at the HE-LHC with 15 ab −1 reaches κ ∼ O(0.1) at 2σ level. Finally, the prospects for directly probing CP violation in the top-quark Yukawa interaction could be potentially further improved by even higher production cross-sections and luminosities offered by the proposed 100 TeV FCC-hh collider [37][38][39], as well as through better background mitigation techniques, especially in the case of thj production, and potential phase-space dependent optimization (reweighing) of CP-odd observables in tth production (see e.g. Ref. [14]), all of which we leave for future work.