Two-loop amplitudes for t ¯ tH production: the quark-initiated N f -part

: We present numerical results for the two-loop virtual amplitude entering the NNLO corrections to Higgs boson production in association with a top quark pair at the LHC, focusing, as a proof of concept of our method, on the part of the quark-initiated channel containing loops of massless or massive quarks. Results for the UV renormalised and IR subtracted two-loop amplitude for each colour structure are given at selected phase-space points and visualised in terms of surfaces as a function of two-dimensional slices of the full phase space.


Introduction
Higgs production in association with a top quark pair was observed for the first time a few years ago at the Large Hadron Collider (LHC) [1][2][3] and will play an important role at the High-Luminosity (HL) LHC.The process pp → t tH is particularly interesting due to its direct sensitivity to the top-Yukawa coupling y t , which is now being constrained with increasing accuracy, including potential CP-violating couplings [4,5].The importance of this process was realised a long time ago [6,7], and NLO QCD corrections for on-shell t tH production have been known for many years [8][9][10][11][12].The corrections have been matched to parton showers in Refs.[13][14][15].NLO EW corrections have first been calculated in Ref. [16], the EW corrections have been combined with NLO QCD corrections within the narrow-width-approximation (NWA) for top-quark decays in Refs.[17,18].NLO QCD corrections to off-shell top quarks in t tH production with leptonic W -decays have been calculated in Ref. [19,20] and full off-shell effects with H → b b have been calculated in Refs.[21,22].A combination of the NLO QCD corrections with NLO EW corrections has been presented in Ref. [23], NLO QCD corrections combined with electroweak Sudakov logarithms and a parton shower have been studied in Ref. [24].
Given the projection that the statistical uncertainty will shrink to the order of 2-3% after 3000 fb −1 [33], the measurement of t tH will be dominated by systematics.As the dominant systematic uncertainties currently come from modelling uncertainties of signal and backgrounds [1][2][3], there is a clear need to reduce the theory uncertainties.At NLO QCD the scale uncertainties are of the order of 10-15%, therefore NNLO QCD corrections are necessary to match the experimental precision at the HL-LHC.
First steps towards this goal already are available in the literature: in Ref. [34], O(α 4 s ) corrections to the flavour-off-diagonal channels have been calculated, exploiting relations from q T -resummation [35].In Ref. [36], the total NNLO cross section has been presented, where for the finite part of the two-loop virtual amplitude a soft Higgs boson approximation has been used.The coefficients of the two-loop infrared singularities for this process have been calculated in Ref. [37].In Ref. [38], the order O(y 2 t α s ) corrections to the perturbative fragmentation functions and to the splitting functions relevant for associated top-Higgs production have been calculated.Analytic results for the master integrals entering the leading-colour two-loop amplitudes that are proportional to the number of light flavours for the processes gg, q q → t tH have recently been presented in Ref. [39].Furthermore, the gg → t tH one-loop amplitude has been calculated semi-numerically up to second order in the ε-expansion in Ref. [40].Results for the two-loop amplitudes for both the gluon and the quark channel in the high-energy boosted limit have been provided very recently in Ref. [41].
In this work we present numerical results for the two-loop virtual amplitudes for q q → t tH which contain closed fermion loops, i.e. are proportional to the number of light fermion flavours n l , heavy fermion flavours n h , or both.Specifically, we calculate the renormalised interference of the two-loop amplitude with the tree-level amplitude, with full dependence on the top quark and Higgs masses, split into nine independent colour and fermion flavour factors.Many of the master integrals appearing in this calculation are not currently known fully analytically, we therefore choose to evaluate all integrals using the sector decomposition [42][43][44][45] approach.Our results are visualised on one-and twodimensional slices of the five-dimensional phase space.These results can be regarded as a proof of concept for the calculation of the other colour structures and partonic channels.
The paper is structured as follows.In Section 2 we describe the kinematics of the process, the structure of the q q → t tH amplitude, and the workflow of our calculation.We present our results in Section 3 and conclude in Section 4. Further details of our calculation, including the UV renormalisation, the colour decomposition, the integral families used for the integral reduction, and full numeric results at several example phase-space points are given in Appendix A.

Description of the method
The calculation of the virtual two-loop amplitudes contains the channels q q → t tH and gg → t tH.Here we focus on the quark initial state.

Kinematics
We use "all incoming" kinematics, and the Mandelstam invariants are defined by Ten such invariants can be built; five of them are independent due to momentum conservation.Out of these we use the following dimensionless variables: There is also an independent parity odd invariant, The square of the parity odd invariant is equal to the Gram determinant spanned by four linearly independent external momenta and is not algebraically independent of the other invariants.However, as ϵ(1234) picks up a sign under parity, while the square root of the Gram determinant does not, the sign of the parity odd invariant must be specified to fully describe a physical phase-space point.QCD is invariant under parity, therefore, the QCD corrections to the t tH production amplitudes ultimately must not depend on the sign of the invariant.
The phase-space parameters.The angles θ t and φ t are local to the t t rest frame, while θ H is local to the t tH rest frame.

Phase space parametrisation
The phase-space volume for t tH production is non-trivial when expressed in the variables given in eq.(2.4).To parametrise it in a more explicit way we factorise it into sub-phasespace volumes for the production of a "t t state" and a Higgs boson, combined with the "decay" of the t t state into two top quarks, leading to the following expression: and the angles θ H , θ t , and φ t introduced as in Figure 1 with precise definitions given in Appendix A.4.
As the production threshold of the t tH system is located at s 0 = (2m t + m H ) 2 , a convenient variable for a scan in partonic energy is such that β 2 → 0 at the production threshold and β 2 → 1 in the high energy limit.For a compact parametrisation of the fraction of kinetic energy which enters the t t system, we define the variable frac s t t as with frac s t t = 0 corresponding to the production threshold of the t t system with the Higgs boson carrying the remaining energy, and frac s t t = 1 corresponding to the production threshold of the Higgs boson, with the t t system carrying the remaining energy.Note that if the phase-space integration is performed in frac s t t , a Jacobian factor of ds t t/dfrac s t t has to be included for the full phase-space density of (2.10) The set of parameters {β 2 , frac s t t , θ H , θ t , φ t } provide a way to parametrise the amplitude which is equivalent to using the five invariants from eq. (2.4); the mapping between them is defined by eq.(2.7), eq.(2.8), eq.(2.9), and the relations given in Appendix A. 4. In these parameters the physical region of the phase space is found as (2.11) Note that the probability density of eq.(2.10) will suppress the low-β region as β 4 , and enhance the high-β region as 1/(1 − β) 2 .It will also suppress both low-and high-frac s t t regions as frac s t t and 1 − frac s t t , respectively.Nominally, the factors sinθ H and sinθ t also suppress the polar cap regions in θ H and θ t , but this is only an artifact of the choice to map the respective spherical regions to a hypercube.

Ultraviolet renormalisation
To produce an ultraviolet (UV) and infrared (IR) finite 2-loop amplitude, we work in conventional dimensional regularisation assuming that all momenta live in d = 4 − 2ε space-time dimensions, and use the expressions for the two-loop singularity structure of massive amplitudes worked out in Ref. [46], also used in Ref. [47].
We expand the bare amplitude for q q → t tH in the strong coupling as where the dependence on kinematics is implicit.The renormalised amplitude has the form where Z q (Z Q ) are the on-shell wave-function renormalisation constants for light (heavy) quarks, µ is the renormalisation scale, and the bare quantities are replaced by the respective renormalised ones.The bare mass of the heavy quark, m 0 , is renormalised using Z m in the on-shell scheme and the bare Yukawa coupling y 0 t is renormalised accordingly.In particular, we further expand the coefficients of the expansion of the bare amplitude in eq.(2.12) to where and the overall factor of m −1−nε is extracted in order to have dimensionless amplitudes that depend on the dimensionless variables x ij introduced in eq.(2.4).
The bare coupling constant α 0 s is defined as which corresponds to the MS scheme with N f = n l + n h active flavours.However, we do not consider the top quark as an active flavour contributing to the running of α s and the parton distribution functions.Therefore we use the decoupling relation where ζ αs is given in Appendix A.1, together with the explicit expressions for the renormalisation constants.Our notation is such that α s ≡ α . To split two-loop amplitudes into smaller building blocks, it is in general convenient to project the amplitude onto scalar form factors multiplying the independent spinor and Lorentz structures, or onto helicity amplitudes.The latter is however less convenient for amplitudes involving massive fermions.For the q q channel considered here, we use the Born amplitude itself as a projector, and calculate the spin-and colour-summed interference term of the renormalised and rescaled NNLO amplitude with the LO amplitude.We decompose this quantity into colour and flavour factors as follows: where N C is the number of QCD colours, and the colour group factors C F , C A , T F , and d 33 are given in Ref. [48] for SU (N ) as well as SO(N ) and Sp(N ): we allow the colour group to be general here.
Similarly, the decomposition of the NLO and LO interference terms is ) The explicit expressions for the renormalized components A, B i and C i in terms of their bare counterparts are given in eq.(A.9) to eq. (A.12) in Appendix A.1.Note that the colour factors of the t tH production amplitude are in principle the same as for t t production.The virtual amplitudes for top quark pair production at NNLO were investigated in Ref. [47], see also Refs.[49][50][51] for analytic results in the quark channel.In Ref. [47], the colour factor decomposition is given after having formed the interference with the Born amplitude, just as in eq.(2.19).In contrast to [47, eq. (2.8)] however, we do not assume the colour group to be SU (N ), and as a result, instead of seven independent colour and flavour structures, we identify nine.

Infrared singularity structure of the virtual amplitude
The origin of the IR divergences present in 2-loop amplitudes with two massive coloured final state particles has been discussed in the literature and their form at 2-loops is known [46,47,52].
For the description of the IR divergences, we work in the colour space formalism.The renormalised amplitude is expressed as a vector in colour space |A R (α s , y t , m, µ, ε)⟩ and the divergences are removed by using a multiplicative MS renormalisation factor Z ≡ Z({p}, {m}, µ, ε), which is a kinematic-dependent matrix in colour space.For the colour decomposition for q l qk → t i tj H we adopt the following basis elements1 : with full details given in Appendix A.2.We use the following expressions: The renormalisation constant Z fulfills the differential equation which induces the following solution at two loops: 4ε . (2.27) The coefficients of the anomalous dimension matrix are defined by the expansion with The general form of the anomalous dimension matrix up to two loops is given in Ref. [46].Here we present the explicit form of the expressions for the process we study.With the colour basis of eq. ( 2.24) the anomalous dimension matrix has the specific form with and in terms of light flavour and colour factors leads to 11, n l d 33 + . . .
where the parts irrelevant for the N f -part of the amplitude are contained in the dots.The relevant Z 11 components then have the form which can now be used in order to construct the IR structure of the interference terms, as will be shown in the following: The parts not shown here have no IR poles.

Workflow of the calculation
The leading order (LO) amplitude A b 0 can be represented by two Feynman diagrams: (2.36) The LO amplitude has no N f -part itself, but it contributes to the renormalisation of the NNLO N f -part, because the α s beta-function contains N f .The LO amplitude in the quark channel has both ε 0 and ε 1 parts (but no higher parts).We derive the corresponding expression using Alibrary [53], which is a Mathematica library interfacing with Qgraf [54], Feynson [55, Chapter 4], Form [56], and Color.h[48] to generate amplitudes, sum over tensor structures, construct integral families, and export the results to integration-byparts (IBP) relation solvers and/or pySecDec [57][58][59][60].
We can use the LO result to estimate the distribution of the events over the phase space at the LHC, as done in Figure 2.These plots tell us that most of the events are expected to come from the region of moderately high β2 and medium frac s t t .In particular, the region of β 2 ∈ [0.10, 0.95] (that is, √ ŝ ∈ [500 GeV, 2.1 TeV]) is expected to contain 99% of all events.

Amplitude generation
To generate the one-loop and two-loop amplitudes (A b 1 and A b 2 respectively) we use the following procedure: first we generate the corresponding Feynman diagrams (using Qgraf), then we insert Feynman rules, apply the projectors, and sum over the spinor and colour tensors (using Form and Color.h);all of this is done through Alibrary.This way, for each diagram, we obtain a corresponding sum of many scalar integrals.
In total we find 31 non-zero one-loop diagrams and 249 two-loop diagrams.Examples for one-loop diagrams with different colour factors are depicted in Figure 3; examples for two-loop diagrams can be found in Figure 4. 2

IBP reduction
The next step is to reduce the calculation of the approximately 20000 scalar integrals that appear in the amplitudes to a much smaller number of master integrals using IBP relations [65].To this end we first calculate the symmetries between the diagrams (using (left), and β 2 and frac s t t (right), according to the LO q q → t tH amplitude.For this plot we take the energy of incoming quarks to be distributed according to the ABMP16 parton distribution functions [61] (which we evaluate via LHAPDF [62]), with the collision energy set to 13.6 TeV.We have also applied cuts on the top quark momenta (as we calculate with on-shell top-quarks) in line with those reported in [1,3]: we enforce a minimal transverse momentum of 25 GeV, a maximal rapidity of 4.5, and a separation ∆R in rapidity and azimuthal angle between the top quarks of ∆R > 0.4.These cuts remove about 3% of the events, and mostly affect the low-β region.
Feynson), and sort them into integral families.For the q q channel we use 4 one-loop and 43 two-loop families.Out of the two-loop families, 28 are unique (shown in Appendix A.5); the remaining 15 are obtained by crossing the external momenta.
The usual next step would be to write down a system of IBP equations and solve it symbolically using e.g. the Laporta algorithm [66].We employed Kira [67,68] together with Firefly [69,70] for this purpose and observed that while the reductions for the oneloop families may be obtained in a rather straightforward manner, for most of the two-loop families the computation is quite challenging.Fully analytic reduction of the two-loop amplitude is rendered largely unfeasible given the large number of variables (five ratios given in eq.(2.4), the mass ratio m 2 H /m 2 t , and the dimensional regulator ε) along with the presence of internal masses.Instead we opt for a numerical approach and solve the IBP systems for individual phase-space points by substituting kinematic scales with rational numbers.Note that we employ the same numerical approach for one-loop amplitudes as well so as to have a uniform implementation for the whole calculation.
To set up the IBP reduction, we first select a basis of master integrals for each of the amplitudes.We require a total of 33 master integrals for the one-loop amplitudes and 831 master integrals for the two-loop amplitudes.The choice of master integrals significantly impacts both amplitude reduction as well as numerical evaluation of the master integrals.Ideally we prefer a basis that is 1) quasi-finite [71], 2) d-factorising [72,73], 3) fast Example diagrams for q q → t tH at one-loop level.Massive quarks are depicted using solid (blue) bold lines, while massless quarks are represented by lighter (grey/red) solid lines.The colour factors correspond to applying the first colour projector from eq. (2.24).
to evaluate with pySecDec, and 4) leads to simple polynomials in the denominators of the IBP reduction tables.Finding a basis satisfying the first two requirements is rather straightforward by considering integrals in higher dimensions (d = 6 − 2ε or d = 8 − 2ε) and with higher propagator powers or dots, with up to 5 dots in some cases.We then apply heuristic arguments to choose integrals that also satisfy the last two requirements in the following way.For each sector, we perform a reduction neglecting sub-sector integrals, and we analyze the denominator factors of the resulting IBP tables for different choices of master integrals.Selecting the master integral basis with the smallest denominator factors, we observe a significant improvement in the run-time for the full reduction.With this initial basis choice, we then evaluate the amplitude as discussed below and we identify the integrals with a significant contribution to the evaluation time.The basis of the corresponding sectors is then further refined by repeating the above procedure, restricting the set of candidate masters to integrals showing a relatively fast convergence with pySecDec.
After we select the basis, we use Kira to generate the IBP systems for each integral family.We generate dimensional recurrence relations using Alibrary to be able to reduce the amplitudes to master integrals in shifted dimensions.The combined system of equations is then fed to Ratracer [74] which prepares and optimizes an execution trace of the solution.Then we use Ratracer to perform a series expansion on this trace in ε; this results in direct output of the ε-expansion of the IBP coefficients.Performing an expansion in ε effectively removes it from the computation which, combined with substituting rational numbers for the kinematic scales, results in a purely numerical system.This system is then Example diagrams for q q → t tH at two-loop level proportional to n l or n h .Massive quarks are depicted using solid (blue) bold lines, while massless quarks are represented by lighter (grey/red) solid lines.The colour factors correspond to applying the first colour projector from eq. (2.24).
solved by Ratracer through replaying the trace in a parallelized manner and using finite field methods.Note that finite field methods used for function reconstruction as a way of solving IBP equations is by now an established practice, pioneered in Refs.[75,76]; our usage however does not require function reconstruction, only rational number reconstruction and the Chinese remainder theorem.Our setup allows us to compute reductions in around two CPU minutes for the two-loop amplitude, and under a second for the one-loop amplitude on a desktop CPU for most points.Overall this reduction method is fast enough, in the sense that we are more constrained by the evaluation of the master integrals.

Evaluation of the master integrals
The families of integrals required in this calculation are complicated enough that analytic evaluation is not currently achievable for many of the master integrals, such as topologies b81 or b82 shown in Appendix A.5.Instead, we turn to evaluating the master integrals numerically, using the approach of sector decomposition as implemented in pySecDec.Specifically, we rely on pySecDec's ability to integrate weighted sums of integrals (introduced in version 1.5, see [58]), and define one sum for each of the bare amplitude's symbolic structures as given in eq.(2.19) and eq.(2.21).We require the two-loop amplitudes to be expanded up to ε 0 , and one-loop up to ε 1 .By default pySecDec expects the coefficient of the sums to be given as algebraic expressions in terms of kinematic variables, but because we do not compute these in general (as we do not perform a fully symbolic IBP reduction), we instead use the results of a test IBP reduction at some fixed kinematic point for the coefficients to compile the pySecDec integration library.This ensures that the pole structure of the amplitude is known at the compilation stage, and so the needed depths of the ε-expansion of the masters can be correctly determined.After the integration libraries are compiled, to evaluate the amplitudes at a given kinematic point, we substitute the coefficients with the results of the IBP reduction.
The sector decomposition of the 831 two-loop master integrals results in a total of around 18000 sectors, and around 28000 sector/expansion-order pairs.To make the evaluation of such a lage number of expressions efficient we rely on the performance improvements in pySecDec 1.6 (see [57]) coming from the new Disteval evaluator (which was partially developed as a response to the challenges of this calculation).We perform all the evaluations of the two-loop amplitudes on NVidia A100 GPUs.
The one-loop amplitudes on the other hand are much simpler (180 sectors in total) and quicker to evaluate; for them we only use CPUs.
Our target is to obtain the renormalised two-loop amplitudes with a precision better than 1%.In the bulk of the phase space this is easily achievable, and the two-loop integration time per point is around 5 minutes; for the one-loop amplitude it is around 10 seconds using 4 CPU threads.This however changes in the high-β region and in the regions near the boundaries of frac s t t and θ t : there we observe large numerical cancellations, both within and between the integrals, that require the evaluation of the master integrals to higher precisions to meet the amplitude precision goals.
These cancellations cause three separate problems: 1.They drive the evaluation times upward, and in principle we expect this growth to be unbounded as β tends to 1 (i.e.ŝ → ∞).
This problem could be mitigated by an asymptotic expansion in large ŝ for the highβ region.However, we will not follow this strategy here, in favour of using a single method for all regions.
2. They require increasingly large Quasi-Monte-Carlo lattice sizes, to an extent where we run into the limitations of precomputed lattices available in pySecDec: the largest such lattice has the size of around 7•10 10 , but some of the integrals need up to 10 14 evaluations in the high-β region.
To solve this issue we employ the new "median QMC rules" [77] lattice construction method implemented in pySecDec 1.6, that enables on-the-fly construction of lattices of unbounded size.
3. At very high β 2 (e.g.0.99) the cancellations between some of the integrals become as large as 20 decimal digits, which means that even evaluating the integrals to the full precision of the double-precision floating point numbers (which is 16 digits) would be insufficient to get any precision for the amplitudes.
Here we find that the integrals that cancel between each other and need high precision are mostly simpler integrals with up to five denominators, most significantly the ones of "sunset" and "ice-cone" [78] type, in various mass and external momenta configurations.Such integrals converge relatively quickly, and obtaining them with more than 20 digits of precision using sector decomposition would be well within our time budget if it was not for the double-precision floating point limitation.
As a solution we have upgraded pySecDec with the ability to dynamically switch from double floating precision to "double-double" for integrals that need it, allowing for the maximum of 32 decimal digits of precision.Our implementation of the double-double arithmetic is based on the methods described in [79,80].We choose this approach instead of the more commonly used quadruple precision floating point numbers (float128) because NVidia GPU compilers do not come with the support for either of them, and our benchmarks show that double-double performs around 2.5 times faster than float128 on a CPU, while being simple enough to be implemented on a GPU.Still, the performance penalty of double-double integration is as high as a factor of 20 on the GPU compared to double.
To cross-check our double-double precision implementation we have also evaluated the sunset and ice-cone type integrals using the series solution of differential equations as implemented in the DiffExp package [81,82] with boundary conditions obtained using the Auxiliary Mass Flow (AMFlow) method [83].We find agreement between our results up to the error reported by our double-double implementation.
Once the integrals are evaluated, the last step is to combine the values of the bare Born, one-loop and two-loop results to values for the renormalised virtual two-loop amplitude as described in Section 2.4, propagating the numerical uncertainties.

Checks
To double-check our calculation we have independently computed the LO and NLO amplitudes via GoSam [84,85] at a number of points, verifying agreement within the reported accuracy.Note that the comparison of the NLO amplitudes requires extra care, because GoSam produces the results in the 't Hooft-Veltman scheme [86], and these need to be converted to get agreement with conventional dimensional regularisation that we use.Regularisation scheme independence of the NLO virtual contribution is only obtained after full IR subtraction [87].In the particular case of interest, the scheme difference of the subtraction term can be traced to the O(ε 1 ) part of the Born contribution which vanishes for the 't Hooft-Veltman scheme.Hence, we can convert the GoSam result to conventional dimensional regularisation by (2.37) The factor 1 − π 2 12 ε 2 is necessary because the convention for MS in GoSam uses (4π) ε instead of S ε of eq.(2.14) as a prefactor.
We have also double-checked the IR poles of our amplitudes against Ref. [37], where the pole parts of the renormalised interference terms are given at four phase-space points.To get agreement with this paper we need to set the renormalisation scale to m t , and use 6 fermion flavours in the running of α s .
While some of the symmetries listed in Section 2.3 are trivially observed when deriving the variables of eq. ( 2.4) from the parameters of eq.(2.11), we have verified the symmetry of simultaneous exchange of q ↔ q and t ↔ t.
Finally and most importantly, for each evaluated point we compute the predicted IR pole coefficients of the amplitude as described in Section 2.5, and compare them to the ones obtained by numerical integration.We find agreement within the reported integration accuracy, which provides us with a check on both the correctness of the renormalisation procedure, and on the correctness of the reported numerical integration precision of the two-loop results, since the IR poles only depend on the one-loop amplitudes, which are integrated separately (and to a higher precision).

Results
In this section we visualise the two-loop amplitude as a function of slices of phase-space variables.To this end we choose the following kinematic point to centre our slices on: β 2 = 0.8, frac s t t = 0.7, cos θ H = 0.8, cos θ t = 0.9, cos φ t = 0.7; ( we also set m 2 H /m 2 t = 12/23, µ 2 = s 12 /4, and m 2 t = 1.The values of the amplitude at this kinematic point is given in Appendix A.6 (along with two other points), where we list both the bare and the renormalised values of each component of the LO, one-and two-loop amplitude (as defined in eq.(2.20), eq.(2.22), and eq.(2.23)).For brevity however we prefer not to plot individual components, but rather the combined C and B values as defined in eq.(2.19) and eq.(2.21).For this we set n l = 5, n h = 1, and the colour group to SU (3), i.e.C F = 4/3, C A = 3, and d 33 = 5/6.
In Figure 5 we have plotted one-dimensional slices of both C and B in β 2 and frac s t t .The plots illustrate the difference in behaviour of the one-and two-loop amplitudes across the parameter space; the two-loop amplitude changes more rapidly, and is on average more negative.
Both in Figure 5 and further it is convenient to use the LO amplitude A as a reference; a slice of A in β 2 and frac s t t is presented in Figure 6.Note that once the phase-space  A × (phase-space density) × 10 3 On the right the amplitude is multiplied by the phase-space density of eq.(2.10).The centre point is marked with a star.density factor of eq.(2.10) is included to obtain the event probability density, the regions of low-β 2 , low-frac s t t , and high-frac s t t are all suppressed.This suppression is important because starting at the one-loop level the amplitude develops a Coulomb-type singularity in the low-frac s t t region.This singularity can be seen on the slices of B and C in β 2 and frac s t t depicted in Figure 7.The inclusion of the phase-space density however suppresses 0.2 0.4 0.6 0.8 C × (phase-space density) × 10 3 this divergence, as can be observed in Figure 8.
To further illustrate the difference in behaviour between the one-and the two-loop level results, we present the slice in θ H and θ t in Figure 9.A similar slice in θ t and φ t is presented in Figure 10.Finally, we illustrate the difference in behaviour between different components of B and C in Figure 11, with a slice in β 2 and frac s t t for each of the individual components, aside from B l , C ll , which are not plotted because their ratio to A is constant.

Conclusions
We have presented numerical results entering t tH production at NNLO QCD, for the quark initiated N f -parts of the two-loop amplitude including loops of both massless and massive quarks.This calculation serves as a proof of concept that our setup is capable of calculating two-loop pentagon amplitudes with internal massive propagators and three massive particles in the final state.We have performed the UV renormalisation and subtraction of IR poles, presenting the finite part of the two-loop amplitude, split into nine different colour structures for a general colour group.
For the reduction to master integrals, we do not attempt to obtain a fully symbolic reduction and instead perform a numerical reduction for each phase-space point leaving the dimensional regulator symbolic.The master integrals are evaluated with a recent version of pySecDec, which has been further extended to support integration over double-double precision integrands, this allows us to obtain stable results also in the high-energy and collinear limits where many digits of the master integrals cancel.The evaluation times vary substantially over the phase space, being of the order of five minutes in the bulk of the phase space, increasing substantially when approaching the β → 1 limit.We do not expect the full quark channel to present further major obstacles within our calculational framework.
Although we have demonstrated that the amplitude can be evaluated with sufficient precision at individual phase-space points, the largest remaining challenge for producing realistic phenomenological applications is to sufficiently densely sample the full 5-dimensional phase-space.One possible way of addressing this obstacle is to supplement the evaluated phase-space points with a reliable interpolation framework that allows data points at any 5-dimensional phase-space point to be provided with sufficient accuracy.This is a challenge for kinematic regions where the amplitude has a very steep gradient, for example in the high-energy region with quasi-collinear configurations.While an interpolation covering the whole phase space is feasible, assessing the associated uncertainties is challenging; this is work in progress.

A.5 Integral families
Generic integral families needed for the NNLO q q → t tH: In this section we provide results for both the renormalised and bare amplitudes at three example points.The first point is a rationalised version of the centre point from eq. (3.1); it is given by

Figure 2 :
Figure2: Event probability distribution in β 2 (left), and β 2 and frac s t t (right), according to the LO q q → t tH amplitude.For this plot we take the energy of incoming quarks to be distributed according to the ABMP16 parton distribution functions[61] (which we evaluate via LHAPDF[62]), with the collision energy set to 13.6 TeV.We have also applied cuts on the top quark momenta (as we calculate with on-shell top-quarks) in line with those reported in[1,3]: we enforce a minimal transverse momentum of 25 GeV, a maximal rapidity of 4.5, and a separation ∆R in rapidity and azimuthal angle between the top quarks of ∆R > 0.4.These cuts remove about 3% of the events, and mostly affect the low-β region.

Figure 5 :
Figure5: One-dimensional slices in β 2 and frac s t t of the one-and two-loop amplitudes B and C as defined in eq.(2.22) and eq.(2.20), normalised to the Born amplitude squared A from eq. (2.23), around the centre point of eq.(3.1).The centre point is marked with a star.Each plot is an interpolation from around 30 data points.

Figure 6 :
Figure 6: Slice of the LO amplitude around the centre point of eq.(3.1) in β 2 and frac s t t .On the right the amplitude is multiplied by the phase-space density of eq.(2.10).The centre point is marked with a star.

Figure 8 :
Figure8: Slices of the one-loop (left) and two-loop (right) virtual amplitudes multiplied by the phase-space density of eq.(2.10), around the centre point of eq.(3.1) in β 2 and frac s t t .The centre point is marked with a star.

A. 6
Numerical results at example phase-space points α 2 = sp 13 sp 14 sp 23 sp 24 m 4 µ 4 , α 3 = sp 13 sp 24 sp 23 sp 14 , β 34 = acosh −sp 34 2m 2 , (2.30) and sp ij ≡ 2p i •p j +i0 + .The anomalous dimensions γ i are given in Appendix A.3.Since we are only interested in the interference with the LO amplitude, we only need the component Γ 11of the anomalous dimension matrix Γ for the IR pole structure of those amplitude parts proportional to the quark flavours at NNLO.Expanding in α (n l ) s frac s t t of the one-and two-loop amplitudes B and C as defined in eq.(2.22) and eq.(2.20), normalised to the Born amplitude squared A from eq. (2.23), around the centre point of eq.(3.1).The centre point is marked with a star.Each plot is an interpolation from around 30 data points.
2Slices of the normalised one-loop (left) and two-loop (right) virtual amplitudes around the centre point of eq.(3.1) in β 2 and frac s t t .The centre point is marked with a star.Each plot is a linear interpolation of grid of around 500 data points in total.
Slice of the normalised one-loop (left) and two-loop (right) virtual amplitudes around the centre point of eq.(3.1) in θ H and θ t .The centre point is marked with a star.Figure 10: Slices of the normalised one-loop (left) and two-loop (right) virtual amplitudes around the centre point of eq.(3.1) in θ t and φ t .The centre point is marked with a star.
Contributions from the individual colour factors to the one-and two-loop amplitudes for phase-space slices around the centre point of eq.(3.1) in β 2 and frac s t t .The centre point is marked with a star.

Table 2 :
Results for the bare amplitudes at the example points from Appendix A.6.