NLO QCD Predictions for off-shell $t \bar t$ and $t \bar t H$ Production and Decay at a Linear Collider

We present predictions for $t \bar t$ and $t \bar t H$ production and decay at future lepton colliders including non-resonant and interference contributions up to next-to-leading order (NLO) in perturbative QCD. The obtained precision predictions are necessary for a future precise determination of the top-quark Yukawa coupling, and allow for top-quark phenomenology in the continuum at an unprecedented level of accuracy. Simulations are performed with the automated NLO Monte-Carlo framework WHIZARD interfaced to the OpenLoops matrix element generator.


Introduction
The top quark is the heaviest particle of the Standard Model (SM), and its detailed study offers great potential to probe the electroweak, flavor and Higgs sector. The close connection between the Higgs boson and the top quark is most apparent for the (meta-) stability of the electroweak vacuum, which crucially depends on m t and m H [1][2][3]. A precise determination of top-quark properties is thus a powerful opportunity to find possible hints of new physics and has far reaching consequences for our understanding of the universe. However, at hadron colliders like the LHC many quantities in the top-quark sector, like the top-quark mass, forward-backward asymmetry or the top Yukawa coupling can only be measured with a limited precision. A future linear lepton collider, such as the proposed International Linear Collider (ILC) [4,5] or Compact Linear Collider (CLIC) [6], on the other hand will reach unprecedented precision in the electroweak and top sector.
With respect to top physics, the two most interesting processes to be studied in lepton collisions are top-pair production with and without an associated Higgs boson. Top-pair production allows to measure the top-quark mass at threshold in a theoretically well defined short distance scheme, like the 1S [7] or PS scheme [8], with uncertainties at or below 100 MeV [9][10][11]. Associated ttH production is our best handle to measure the top Yukawa coupling with per cent level precision, see e.g. Ref. [12,13]. Obviously, these physical parameters can only be extracted with this level of precision when the theoretical uncertainties at least match their experimental counterparts.
Due to the relatively large top width, which comes exclusively from the decay into a bottom quark and a W boson, top quarks decay before they can form bound states. The produced W boson decays further via hadronic or leptonic channels, whereas the bottom quark hadronizes and can be identified as a tagged jet. Especially in the clean lepton collider environment, the charge of the b-jet can be reconstructed with reasonably high efficiency [14]. A consistent treatment of the associated finite width effects is both a conceptionally as well as computationally nontrivial problem. Within the so-called narrow-width approximation (NWA), top quarks are produced on-shell and decay subsequently according to their (spin correlated) branching ratios. Higher-order QCD predictions for on-shell top-pair production are well-known, the current best predictions being N 3 LO [15] at the inclusive and NNLO at the fully differential level [16]. First NLO electroweak corrections have been obtained in Ref. [17]. For top-pair production in association with a Higgs boson, there are comprehensive studies of NLO QCD corrections available in Ref. [18]. First inclusive combined electroweak and QCD corrections have been computed in Ref. [19], followed by an in-depth study in Ref. [20].
While computationally simple, the NWA has the obvious drawback that various nonresonant background processes are not included. For off-shell tt or ttH production, however, especially single-top resonances can contribute significantly and can hardly be distinguished experimentally from double-resonant contributions [21]. Furthermore, off-shell effects can only be treated approximatively via a Breit-Wigner parameterization, as in Ref. [22]. Non-resonant contributions and finite width effects can be consistently taken into account employing the complex-mass scheme [23], which guarantees gauge invariance at NLO -at the price of increased computational complexity. Such a calculation for the process e + e − → W + W − bb at NLO QCD has first been presented in Ref. [24]. It has recently been reevaluated in Ref. [25], with the aim of extracting the top-quark width via ratios of single-to double-resonant signal regions.
In this paper, we study top-pair and Higgs associated top-pair production and decay including non-resonant contributions, off-shell effects and interferences at NLO. The simulation is done with the multi-purpose event generator Whizard [26,27], which has been ex-tended to perform automated NLO calculations. In this framework, we compare the on-shell processes e + e − → tt and e + e − → ttH with the off-shell processes e + e − → W + W − bb and e + e − → W + W − bbH. At the differential level, the full processes including leptonic decays are considered, i.e. e + e − → µ + ν µ e −ν e bb and e + e − → µ + ν µ e −ν e bbH. To our knowledge, NLO studies of e + e − → W + W − bbH or the complete off-shell processes e + e − → bb4f or e + e − → bb4f H have not been performed previously in the literature. In contrast, at hadron colliders off-shell top-pair production has been studied in Refs. [28][29][30][31][32][33]. Furthermore, employing the resonanceaware method of Ref. [34], the process pp → bb4f has been matched consistently to parton showers, as presented recently in Ref. [35]. For hadron colliders, corresponding NLO QCD corrections to top-quark pair production in association with a Higgs boson [36] or a jet [37] including leptonic decays have also been studied.
While at hadron colliders top-pairs originate from QCD production, at lepton colliders they are produced via electroweak interactions. This implies that a fixed-order computation of the off-shell processes at a lepton collider comprises a considerably larger set of irreducible electroweak background processes. Such processes involve (very) narrow resonances, like e.g. H → bb. In NLO computations, resonances with very small widths can severely hamper the quality of the infrared (IR) subtraction and consequently influence the convergence and quality of the integration. In order to have these resonance effects under control, in Whizard we have implemented an automatized version of the resonance-aware scheme of Ref. [34]. This is also a prerequisite for a future consistent matching of off-shell processes with parton showers.
Besides the phenomenological relevance of the presented results -in particular for top quark mass measurements in the continuum and measurements of the top Yukawa couplingthis paper demonstrates the progress on Whizard as a fully automated NLO event generator. Whizard has for a long time been a (high-multiplicity) tree-level event generator, where besides its usage in all areas of lepton collider physics, its focus on hadron colliders had been mostly on beyond the Standard Model (BSM) physics. NLO QCD corrections have only been considered for the explicit study of pp → bbbb [38,39]. Furthermore, NLO QED effects have been studied on fixed order as well as by resumming soft photons for chargino production at the ILC [40,41]. Apart from Whizard, various collaborations are including generic NLO simulations into their event generators. This has been made possible by tremendous advances in the automation of the computation of one-loop amplitudes during the last decade. Publicly available one-loop providers (OLPs) such as Helac-1Loop [42], OpenLoops [43], GoSam [44], Recola [45,46] or MadLoop [47] can compute arbitrary virtual matrix elements in the SM, though in practice limited by computing power. Complete NLO QCD support has so far been achieved within the frameworks of Helac-NLO [48], Madgraph5_aMC@NLO [49], Sherpa [50] and Herwig7 [51]. Finally, we want to remark on a topic that is both closely related to fixed order NLO predictions and important for the description of tt and ttH. At threshold, these processes actually require the inclusion of bound state effects that can be treated in non-relativistic QCD. Here, the exchange of soft gluons leads to Coulomb singularities that have to be resummed. In order to take this into account, Whizard ships with the Toppik program [7], which can be used to compute resummed form factors up to next-to-leading-logarithmic (NLL) accuracy, as presented in [52]. Inclusive NNLO calculations at threshold have been compared in Ref. [53]. The state-of-the-art of fixed-order corrections has been recently improved to N 3 LO [54] and at the resummed level to NNLL [55]. However, these results are only accurate in the threshold region, while Whizard can describe both the threshold and continuum domain by using a smooth matching approach. Whizard's NLO capabilities are used hereby to compute the radiative corrections to the top decay in a factorized approach as well as to obtain the full W + W − bb process at NLO. A preliminary status thereof is presented in Ref. [56,57], while an in-depth study of the threshold matching in Whizard is in preparation [58]. In this paper, we focus on a fixed-order description of the continuum, while pointing out regions where threshold effects become important.
The paper is organized as follows. In section 2, we describe the setup of the calculation. In section 3, we address the issue of resonance-aware subtraction and its implementation in Whizard. In section 4, the phenomenology of tt and ttH is briefly reviewed. The employed input parameters, scale choices and phase-space cuts as well as an overview of the performed validations can be found in section 5. The main phenomenological results of this paper can be found in section 6 and 7, where in section 6 we focus on results at the inclusive level, while in section 7 corresponding differential predictions are investigated. We discuss scale variations for the NLO QCD corrections, show results for polarized lepton beams and discuss the influence of the NLO QCD corrections on the extraction of the top Yukawa coupling. Additional differential results as well as technical details on the resonance-aware subtraction can be found in the appendix. Our conclusions are presented in section 8.

Setup of the calculation
The predictions presented in this paper are obtained with the automated Monte Carlo framework Whizard combined with the amplitude generator OpenLoops. As in this paper we introduce this framework for the first time, in the following, we give a short introduction to both programs, starting with a discussion of the event generator Whizard and its treatment of next-to-leading order QCD corrections in section 2.1. This is followed by a description of OpenLoops in section 2.2.

The
Whizard event generator at next-to-leading order Whizard [26,27] is a multi-purpose event generator for both lepton and hadron colliders. At leading-order, it can deal with arbitrary SM processes, as well as a multitude of BSM processes (e.g. generated from automated tools like in Ref. [59]). Moreover, it can perform simulations for a broad class of processes at next-to-leading order. The modern release series (v2) has been developed to meet the demands of LHC physics analysis, while its generic treatment of beam-spectra and initial-state photon radiation makes it especially well suited for lepton collider physics. The program has a modular structure and consists of several subcomponents, the most important being O'Mega [27], Vamp [60] and Circe [61]: O'Mega computes multi-leg treelevel matrix elements as helicity amplitudes in a recursive way that avoids Feynman diagrams. Vamp is used for Monte-Carlo integration and grid sampling. It combines the multi-channel approach [62] with the classic Vegas algorithm [63] to automatically integrate cross sections with non-factorizable singularities. The Circe package can be used to create and evaluate lepton beam spectra.
Whizard can be used for event generation on parton level as well as for the subsequent shower and hadronization. For this purpose, it has its own analytical [64] as well as k T -ordered parton shower, and a built-in interface to Pythia6. Color information is treated in Whizard using the color-flow formalism [65].
The generic NLO framework in Whizard builds upon the FKS subtraction scheme [66], which partitions the phase space into regions where only one divergent configuration is present. This divergence is then regulated using plus-distributions. FKS subtraction allows for the application of Whizard's optimized multi-channel phase-space generator for the underlying Born kinematics, from which real kinematics are generated. It is also very well suited to the matching procedures employed, as described below. First preliminary results of the NLO functionality of Whizard have been presented in Refs. [57,67]. Whizard supports OpenLoops and GoSam (an interface to Recola is being developed) as one-loop matrix element providers as well as for the computation of color-and spin-correlated Born matrix elements. At tree-level, they can also be used as alternatives to O'Mega.
For event generation, Whizard can produce weighted fixed-order NLO QCD events that are written to HepMC [68] files. This allows for flexible phenomenological fixed order studies, especially in combination with Rivet's [69] generic event analysis capabilities. Matching to parton showers is achieved with an independent implementation [70] of the Powheg matching method [71].
Apart from scattering processes, Whizard is also able to compute decay widths for 1 → N processes at NLO. The final-state phase space is built in the usual fashion, whereas the initialstate phase space is adapted for decays. Whizard constructs the gluon momentum separately and then applies a recursive reassignment of the virtualities of the intermediate particles. Computing decay widths directly in Whizard allows for a consistent treatment of the top width, which has to be recomputed according to the physical parameters and the process definition as discussed in section 5.1.

Virtual matrix elements from OpenLoops
All necessary Born and one-loop amplitudes together with the color and helicity correlators required within the FKS subtraction are provided by the publicly available OpenLoops program [72]. It is based on a fast numerical algorithm for the generation of Born and one-loop scattering amplitudes by means of a hybrid tree-loop recursion that generates cut-open loops as functions of the circulating loop momentum [43]. Combined with the CutTools [73] OPP reduction [74] library and the OneLOop library [75] or with the Collier [76] tensor integral reduction library based on Refs. [77][78][79], the employed recursion permits to achieve very high CPU performance and a high degree of numerical stability. A sophisticated stability system is in place to rescue the small number of potentially unstable phase space points via a re-evaluation at quadruple precision.
Within OpenLoops, ultraviolet (UV) and infrared (IR) divergences are dimensionally regularized and take the form of poles in (4−D). However, all ingredients of the numerical recursion Figure 1: Example pentagon diagrams contributing to the W + W − bb final-state process containing one or two (leftmost diagram) top resonances and an hexagon diagram contributing for W + W − bbH production.
The strong coupling constant is renormalized in the MS scheme, and heavy quark contributions can be decoupled via zero-momentum subtraction in a flexible way, depending on the number of active flavors in the evolution of α S . Unstable particles with a finite width are by default treated via an automated implementation of the complex-mass scheme [23].
The publicly available OpenLoops amplitude library includes all relevant matrix elements to compute NLO QCD corrections, including color-and helicity-correlations and real radiation as well as loop-squared amplitudes, for more than a hundred LHC processes. Many libraries for lepton collisions can easily be taken from this LHC library, as any crossing of external particles is automatically done when a library is loaded. For example, the one-loop library to be used for the process e + e − → jj is ppll. For many other processes, especially for those with massive quarks in the final state, dedicated lepton collider libraries have been added to the public OpenLoops amplitude repository, which will be further extended in the near future 1 .
The Whizard+OpenLoops interface is based on the BLHA standard [88]. Moreover, this new interface can use a modification of this standard allowing the computation of polarized amplitudes. To this end, the process registry can contain dedicated entries for each polarization configuration of initial or final state particles. This implements an automated NLO setup which allows to study effects of beam polarization -an important feature at future linear colliders like the ILC. Table 1 lists information about the computational complexity with respect to the one-loop amplitudes of the processes studied in this paper. Note that the total number of diagrams is not decisive for the computational effort in the OpenLoops recursion formalism. Instead, the crucial point is the maximal number of n-point functions involved. For the bb(W → lν)(W → lν) processes discussed in this paper, the most complex integrals stem from pentagon diagrams, examples for which are depicted in fig. 1. Also shown in fig. 1 is a hexagon diagram contributing to the associated Higgs production process also discussed in this paper. Concerning the complexity of the amplitudes, due to the reduced number of contributing helicity structures the calculation of the off-shell processes including leptonic decays are less involved compared to the corresponding processes with on-shell W-bosons -despite the increased number of diagrams. 3 Resonance-aware FKS subtraction The standard approach to compute automated NLO corrections can be very inefficient if QCD radiation off partons originating from the decay of a resonance is present. In our case, this issue arises from resonant subprocesses of type Z/H → bb and t → W b. As discussed for the first time in Ref. [34], the problem is due to the fact that the momentum of the resonant particle can be different in the Born phase space and the corresponding real phase spaces with one additional gluon momentum 2 . The real-subtracted contribution to the NLO matrix element contains N + 1-particle matrix elements with corresponding kinematics, as well as Born matrix elements with factorized kinematics in the subtraction terms. In collinear and/or soft regions, it is crucial that both terms agree well. However, the presence of resonances significantly affects their cancellation, and hence the convergence of the integration.
To understand this more in-depth, consider the H → bb splitting with the very narrow Higgs resonance Γ H = O(1 MeV). This occurs as a Higgsstrahlung background process to e + e − → W + W − bb and its decays. Thus, the squared matrix element of the total process contains a term with the contribution of the squared Higgs propagator, wherep 2 bb denotes the invariant mass of the bb-pair in the Born phase space. The Higgs propagator in the corresponding real squared matrix element takes the form where now the Higgs virtuality is made up by the invariant mass of the bb-system and the additional gluon, p 2 bbg . Let the change of the Higgs virtuality from the Born to the real phase space be described by ∆ bbg , such that The explicit form of ∆ 2 bbg does not depend on the process, but on the subtraction scheme. In the FKS approach, the real phase space is constructed in such a way that the invariant mass of the recoiling system and the emitter-radiation system are conserved separately. Thus, ∆ 2 bbg consists of boosts and projections of the Born momenta.
Either way, we can define ε =p 2 bb − m 2 H and, for the ratio of weights associated with the H → bb resonances in the emission matrix element and in the related subtraction terms, we find For the real and subtraction terms to match, it is required that D ≈ 1 in the soft as well as the collinear limit. At the resonance, ε → 0, we see that this condition is fulfilled if ∆ 4 bbg m 2 H Γ 2 H . We immediately see that this poses a problem in the collinear limit, since ∆ 4 bbg can become large if a hard-collinear gluon is emitted. However, also in the soft limit a significant mismatch can occur if the denominator m 2 H Γ 2 H is sufficiently small. This is definitely the case for H → bb, with m 2 4 . As already noted in section 1, the problem that H → bb is contained in the off-shell tt process is unique to the lepton collider, as here at LO the production is of s ) at hadron colliders. For our study, we have addressed the problem of narrow resonances implementing the modified FKS subtraction procedure presented recently in Ref. [34] for generic processes in Whizard. This implementation is briefly outlined in the following. More in-depth information and validation can be found in appendix A.
In the so-called resonance-aware FKS approach, in addition to being partitioned into distinct singular regions, the phase space is also separated into resonance regions, according to the resonance structures of the process. In each extended singular region, the real phase space is constructed in such a way that the invariant mass of the particles which originate from the same resonance is kept fixed. In this way, the shift ∆ bbg in eq. (4) is exactly zero by construction, and hence D = 1. This approach makes use of modified FKS mappings which are evaluated in the rest frame of the corresponding resonance. This leads to the problem that the sum over all singular regions does not reproduce the full real matrix element any more. As shown in Ref. [34], this can be solved by introducing a new component to the integration, the so-called soft mismatch. In Whizard, the integration of the soft mismatch is automatically performed and is included as an additional contribution next to Born, real and virtual components when the resonance-aware FKS subtraction is activated. Related technical details can be found in Appendix A.
In the resonance-aware FKS approach, the standard FKS projectors S α are extended by resonance projectors P α , with P α → 1 if the phase space is close to the resonance associated with the resonance history α . They thus map out this particular resonance structure. Motivated by the narrow-width limit of a resonant process, P α is proportional to the Breit-Wigner factors of a given resonance structure.
In Whizard, resonance information is generated for every simulation, already at leading order. This information is used by the multi-channel integrator Vamp, where all relevant resonance structures are sampled in order to enhance the performance. We use exactly these resonance structures to set up resonance-aware FKS subtraction. Thus, in principle, each of Whizard's integration channels could be identified with the resonance histories, also using the internal mappings used in the construction of the Born phase-space. However, we decided to introduce resonance histories using the projectors of [34] completely independent of the Monte Carlo integration channels.
The implementation of the resonance-aware FKS subtraction led to a restructuring of the HepMC output of weighted fixed-order NLO events. In earlier Whizard versions [67], different phase space points and weights were assigned to each singular region α r . However, when resonances are included, different α r can be associated with the same real phase space (e.g. in the case of Z/H → bb), which leads to an unnecessary abundance of real-emission events in the event output. Therefore, in the most recent Whizard version, a real-emission event is created for each distinct phase-space structure, which is defined by its emitter and the decaying particles. Each of these phase-space structures can be associated to multiple singular regions, over which it is summed to obtain the complete real weight of the event. Moreover, the soft mismatch is included in the subtraction weight.
Employing the resonance-aware FKS subtraction scheme for off-shell top-pair production and decay in leptonic collisions is not trivial, since resonance histories where the gluon is emitted from the production process, i.e. from the top before it decays, cannot be associated to any valid FKS sector. In proton-proton collisions, the consistent resonance-aware treatment of such resonance structures, which requires mappings that preserve simultaneously the invariant masses of the W + b and W −b pairs without the emitted gluon, is guaranteed through FKS sectors associated with the initial-state quark or gluon emitters. However such FKS sectors are not present in the case of uncolored initial states. Thus the extension of the resonance-aware approach to e + e − collisions requires a dedicated treatment for the case of QCD radiation that is emitted by unstable colored particles before they decay.
While this issue deserves more detailed studies that we have deferred to the future, for the study of the off-shell processes e + e − → W + W − bb(H) and e + e − → µ + ν µ e −ν e bb(H), presented in this paper, we use only the two resonance histories Z → bb and H → bb. Employing the implementation of the resonance-aware subtraction scheme with these resonance histories we observe a decent convergence of the numerical integration at the inclusive and differential level.
Finally we want to note, that the resonance-aware FKS subtraction scheme, including a definite resonance history assignment in the event output, enables a consistent matching of fixed-order NLO predictions with parton shower generators for processes with intermediate resonances [34,35]. To this end, all relevant resonance histories should be taken into account.
4 Phenomenology of tt and ttH production and decay

Phenomenology of tt production and decay
In this study we want to investigate NLO QCD perturbative corrections in top-quark pair production at lepton colliders modeling off-shell and interference effects at increasing levels of precision. To this end we will consider the following related 2 → 2, 2 → 4 and 2 → 6 processes, where we treat the bottom quarks as massive. Top quarks almost exclusively decay via t → bW + , such that the process in eq. (6) can be understood as the top-quark pair production process of eq. (5) including top-quark decays. Beyond the narrow width approximation, i.e. including off-shell effects for the produced top quarks, the process of eq. (6), however, receives besides doubly-resonant (signal) top-quark contributions, also contributions from non-resonant and single-resonant (background) diagrams together with their interference. Example diagrams for all three production mechanisms are shown in fig. 2. The sub-dominant single-top diagrams always occur via a fermion line between the two external bottom quarks. Thanks to the finite bottom mass even non-resonant contributions from diagrams with a γ → bb splitting, like the one in the top right of fig. 2, can be integrated over the whole phase space without the necessity for cuts. At the NLO QCD level, the calculation of the process in eq. (6) includes corrections to topquark pair production and also to the top decays together with non-factorizable corrections, which are formally of the order of O (α S Γ t /m t ). Diagrammatically such non-factorizable contributions interconnect production and decay stage of the signal process or the two individual decays, as for example depicted in fig. 1(left). At the same time NLO interference effects with single-resonant and non-resonant contributions and also spin correlations in the top decay are consistently taken into account.
In order to make contact with experimental signatures and to further increase theoretical precision, the process in eq. (7) introduces -beyond the top-quark decays -also leptonic decays of the W-bosons including respective off-shell effects. Due to the purely EW nature of the leptonic W-boson decays, from a perturbative point of view these additional decays do not increase the computational complexity compared to the process with on-shell W-bosons, i.e. the one of eq. (6). However, besides the more involved phase space integration, the number of contributing diagrams increases substantially due to additional single-and non-resonant contributions, as illustrated in fig. 3. Notabene, in the case of decays with initial-state lepton flavor, diagrams like the one on the right of fig. 3 show a singularity and can not be integrated over the whole phase space without cuts. For brevity, here we focus on the different lepton flavor case but an analysis for the very similar same flavor case can easily be performed with the publicly available Whizard+OpenLoops framework. Hadronic top-quark decays will be investigated in the future. Figure 2: The double-resonant signal diagram (top left) besides example non-resonant (top right) and s-and t-channel single-top diagrams (bottom left and right, respectively) of the process   Figure 6: A representative non-resonant diagram contributing to W + W − bbH production via a quartic ZZHH-coupling.

Phenomenology of ttH production and decay
Similar to top-quark pair production, we consider the following related 2 → 3, 2 → 5 and 2 → 7 processes for the associated production of a Higgs boson together with a top-quark pair with increasing level of precision with respect to off-shell, non-resonant and interference effects, where again all b-quarks are treated as massive.
The diagrams involved in these process are very similar to those of the corresponding tt production processes, apart form the additional Higgs boson that couples now to all massive internal or external particles (t, b, W ± , Z, H). Already on the level of the on-shell processes of eq. (8) this results into two competing contributions, as depicted in fig. 5. The diagram on the left of fig. 5 is proportional to the top Yukawa coupling y t and will be denoted as ttH signal contribution, while the diagram on the right can be considered as irreducible Higgsstrahlung background in the ZH channel with an off-shell Z * → tt decay.
Furthermore, at the level of the off-shell processes of eqs. (9) and (10) -besides topologies already present for the corresponding tt processes with an additional attached Higgs boson -, new contributions arise from quartic EW couplings as illustrated in fig. 6. In such contributions, as before, the tiny width of the Higgs boson requires a resonance-aware subtraction scheme to yield a converging integration at NLO over the whole phase space.

Input parameters, scale choice and phase-space cuts
As input parameters, we use the following gauge-boson, quark and Higgs masses [90], The electroweak couplings are derived from the gauge-boson masses and the Fermi constant, G µ = 1.1663787 × 10 −5 GeV −2 , in the G µ -scheme. The CKM matrix is assumed to be trivial, which is for the most relevant element of our computation (V tb ) consistent with the measured value (1.021 ± 0.032 [90]). Furthermore, using the precisely measured value of G µ automatically absorbs important electroweak corrections to the top decay. For the strong coupling constant we use α s (m Z ) = 0.1185 and a two-loop running including n f = 5 active flavors.
With this setup, the gauge boson and top widths are computed directly with Whizard at LO and NLO. In the NLO computation, we use the mass of the decaying particle as renormalization scale. The obtained LO and NLO gauge boson widths are In our calculation we use Γ Z and Γ W at NLO throughout, i.e. also for off-shell cross sections at LO. This ensures that the effective W and Z leptonic branching ratios that result from e + e − → bb4f (H) matrix elements are always NLO accurate. In contrast, in order to guarantee that t → W b branching ratios remain consistently equal to one at LO and NLO, off-shell matrix elements and the top-decay width need to be evaluated at the same perturbative order. For the top width we employ two distinct sets of values: one for the on-shell decay t → W + b and one for the off-sell decay t → ff b, as also detailed in Ref. [30]. The value used for the off-shell top decay includes decays into three lepton generations and two quark generations. It also involves the W width, for which we use the previously computed NLO value. The numerical values are The Higgs width is set to Γ H = 4.143 MeV.
In the determination of the off-shell top width and in all calculations presented in this paper, intermediate massive particles are treated in the complex-mass scheme [23]. This leads to a gauge invariant treatment of finite width effects as well as perturbative unitarity [91]. On the technical side, it necessitates complex-valued renormalized masses that imply for consistency a complex-valued weak mixing angle For the electromagnetic coupling in the G µ scheme we set which gives α −1 e = 132.16916. For the on-and off-shell tt and ttH processes that we consider in this paper, the renormalization scale µ R is set to for tt processes m t + m H for ttH processes and Our default scale choice corresponds to ξ R = 1 and theoretical uncertainties are probed by scale variations. This is obviously no complete assessment of the theoretical errors, for the LO e.g. we do not obtain an uncertainty band, but it is our best handle on perturbative QCD uncertainties.
Thanks to the finite b quark mass all on-and off-shell tt and ttH processes considered in this paper can in principle be integrated over the whole phase space. However, for processes with final state electrons or positrons a singularity emerges for small photon energy transfers, as depicted on the right-hand side of fig. 3. To avoid this, we apply a mild phase-space cut for these processes For the definition of jets we employ the generalized k T algorithm (ee-genkt in FastJet) [92,93] with R = 0.4 and p = −1. We tag b/b-jets according to their partonic content and denote them as j b and jb. Similarly, in the on-shell processes e + e − → tt and e + e − → ttH, we identify the top quark with the jet containing a top quark. In the discussion of differential cross sections in section 7 we always require at least two b-tagged jets. 3 No further phase-space restrictions are applied.

Validation
To validate the new automated subtraction within Whizard, we have performed various cross checks. All of the following checks have been performed at the per mil level, i.e. differences are all at the few per mil level and within two standard deviations of the MC integration. The NLO top-quark width computed by Whizard has been cross-checked both with the value used in Ref. [30] and the analytical formulae [94,95] e + e − → tt and e + e − → W + W − bb at √ s = 800 GeV Figure 7: Total cross section for on-shell and off-shell tt production as a function of √ s and µ R . In the lower panels of the left plot, we show the K-factor for tt and W + W − bb in green and red, respectively, as well as the ratio of off-shell to on-shell results for LO and NLO in blue and red.
processes, like e + e − → qq and e + e − → tt, have been validated against analytical calculations. For e + e − → W + W − bb, we have performed an in-depth cross check with various other results and generators. The total cross section corresponding to the study of Ref. [25], therein computed with Madgraph5_aMC@NLO [49], has been reproduced. Moreover, we find excellent agreement between Whizard, Sherpa [50] and Munich 4 for the parameter set given in section 5.1. Note, that both Sherpa and Munich use CS subtraction [89,96], while Mad-graph5_aMC@NLO and Whizard use FKS subtraction [66]. The resonance-aware NLO calculation was validated internally, comparing the result with a computation based on the traditional FKS subtraction (see also appendix A). To this end, we used large widths in order to avoid problems with the traditional FKS approach.
6 Numerical predictions for inclusive cross sections

Integrated cross sections and scale variation
We start our discussion of the numerical results with an investigation of the NLO QCD corrections to inclusive top quark pair-production cross sections depending on the center-of-mass energy √ s of the leptonic collisions. In the left plot of fig. 7 we show inclusive LO and NLO cross sections for the on-shell process e + e − → tt and the off-shell process e + e − → W + W − bb together with the corresponding K-factor ratios, defined as Right above the production threshold √ s = 2m t , both LO and NLO cross sections are strongly enhanced, and in the limit √ s → 2m t the NLO corrections to the on-shell process e + e − → tt diverge due to non-relativistic threshold corrections, which manifest themselves as large logarithmic contributions to the virtual one-loop matrix element. Instead, in the off-shell process e + e − → W + W − bb the Coulomb singularity is regularized by the finite top-quark width, and the NLO corrections remain finite. However, threshold corrections introduce a distinct peak in the NLO corrections at √ s = 2m t with a maximum K-factor of about 2.5. Below threshold the cross section drops sharply, but QCD corrections remain significant. Far above threshold the NLO corrections are rather small for both the on-shell and the off-shell processes. For e + e − → tt, the corrections remain positive for all √ s. In fact, for large center-of-mass energies, the effect of the top quark mass becomes negligible and the corrections approach the universal leptonic massless quark pair-production correction factor α s /π. In contrast, the NLO corrections to e + e − → W + W − bb decrease significantly faster for large center-of-mass energies, are at the per cent level for √ s = 1500 GeV, and come close to zero at √ s = 3000 GeV. This corresponds to the fact that the non-resonant irreducible background and interference contributions grow with energy relative to the tt signal contribution, which receives purely positive corrections. Our results suggest that at √ s = 800 GeV, positive corrections to the signal process and negative corrections to the background are of the same order of magnitude and partially cancel each other. This leads to very small NLO QCD corrections. However, at this level the currently unknown and possibly large NLO EW corrections to e + e − → W + W − bb have to be included as well for reliable predictions. Comparing off-shell to on-shell cross sections, we see that they are about equal at threshold, but at √ s = 800 GeV the off-shell prediction is about 20% larger.
In the right panel of fig. 7 we show the variation for √ s = 800 GeV of the e + e − → tt and e + e − → W + W − bb NLO predictions with respect to the renormalization scale µ R in the interval µ R = [1/8, 8] · m t . Within the error band [m t /2, 2m t ] predictions for tt and W + W − bb with fixed top-quark width, Γ t = Γ t (µ R = m t ), vary at the level of a few per cent, however with an opposite slope. To understand this behavior, we show the scale variation of the off-shell process additionally with a scale-dependent width, Γ t (µ R ). With such a consistent setting of the width according to the input parameters, including µ R , scale variations in the off-shell process are very similar to the on-shell one. We note that the scale dependence in the top width is in principle a higher-order effect, such that both approaches are in principle valid to estimate missing higher order effects by means of scale variations. However, in order to properly recover the narrow width limit the parameter settings for the width in the propagator and the decay part of the matrix element have to match, including the scale setting.
Inclusive cross sections for Higgs associated top-pair production are shown in the left panel of fig. 8. Also here we observe an enhancement of the cross sections with a maximum located at around √ s = 800 GeV, i.e. far above the production threshold at 2m t + m H ≈ 471 GeV, where σ incl. ( √ s = 800 GeV) ≈ 2.4 fb. Again, NLO QCD corrections are sizeable due to nonrelativistic Coulomb enhancements close to the production threshold. For the off-shell process e + e − → ttH and e + e − → W + W − bbH at √ s = 800 GeV Figure 8: Total cross section of on-shell and off-shell ttH production subject to √ s and µ R . Extra panels as in fig. 7.
e + e − → W + W − bbH the corrections reach +100% and remain large but finite below threshold, while for the on-shell process they diverge close to threshold. Around the maximum of the cross sections, NLO corrections vanish for both, the on-shell and the off-shell process. Above this maximum, the NLO corrections turn negative, yielding corrections at √ s = 3000 GeV of up to −15% for the on-shell process e + e − → ttH and up to −20% for the off-shell process e + e − → W + W − bbH. Again one should also consider how the off-shell cross sections behave relative to their on-shell counterparts. While at LO the e + e − → W + W − bbH cross section decreases considerably slower with energy compared to the on-shell process e + e − → ttH, at NLO the corrections to the off-shell process are more sizeable and negative with respect to the on-shell case, yielding comparable inclusive cross sections for the on-shell and off-shell process. Still, at 3000 GeV the off-shell inclusive cross section is about 20% smaller then the on-shell one.
In the right panel of fig. 8, we display renormalization scale variations at √ s = 800 GeV for Higgs associated top-pair production. For this center-of-mass energy scale variation uncertainties in e + e − → ttH are negligible (induced by vanishing NLO QCD corrections), while in e + e − → W + W − bbH with the standard choice Γ t = Γ t (µ R = m t ) they amount to several per cent in the considered variation band. Similar to the tt case, we also show scale variations taking consistently into account the scale dependence in the top-quark width. Here, the behavior of the off-shell process is very similar to the on-shell one.
Finally, in tables 2 and 3 we list inclusive cross sections for tt and ttH (both on-and off-shell) processes, respectively, for several representative center-of-mass energies. Listed uncertainties are due to scale variations, where we employ the fixed top-width, Γ t = Γ t (µ R = m t ). In section 7 we will continue our discussion of NLO corrections for top-pair and Higgs associated top-pair

0.79
production at the differential level. There we will focus on √ s = 800 GeV, as here cross sections are largest for ttH production, which should offer the best condition for a precise determination of the top Yukawa coupling, as discussed in the following section. While consider this as a viable running scenario for a precision measurement, one should keep in mind that for other energies the NLO QCD corrections will be larger in general, at least at the inclusive level.

Determination of the top Yukawa coupling
A precise measurement of Higgs associated top-pair production allows for the direct determination of the top-quark Yukawa coupling y t at the per cent level [12,13]. This allows -next to the measurement of the ttZ coupling -for decisive probes of many new physics models, as significant deviations from the Standard Model value y SM t = √ 2m t /v are predicted in many such models, e.g. in generic two Higgs-doublet models, the MSSM or composite Higgs or Little Higgs models. A per cent level measurement of y t is feasible at future high-energy lepton  Figure 9: The e + e − → ttH and e + e − → W + W − bbH LO and NLO cross sections as a function of the top Yukawa coupling modifier ξ t = y t /y SM t , as well as a linear fit used to determine the coefficient κ as described in the text, (21). colliders, as the ttH and W + W − bbH cross sections are quite sensitive to y t . The sensitivity of the ttH processes (on-and off-shell) on y t is commonly expressed in terms of [13,97] ∆y t y t = κ ∆σ σ . (21) In this way, the relative accuracy on the measured cross section can directly be related to a relative accuracy on the top Yukawa coupling. Since the y t -dependence of the cross section is approximately quadratic, κ is close to 0.5. More precisely, parameterizing deviations of the top-Yukawa coupling from its SM value as y t = ξ t · y SM t we can write the total cross section as σ(ξ t ) = ξ 2 t · S + ξ t · I + B, where S and B denote ttH signal 5 and background contributions, respectively, while I stands for interference terms. The y t -sensitivity of ttH cross sections can be determined via a linear fit of σ(y t ), which corresponds to Note that whereas B is strictly positive, we can make no statement about the sign of I. Eq. (22) -making the quite general assumption that the signal dominates over the interference, −I < 2S -shows that κ < 0.5 can only be realized via sufficiently large and negative interference contributions, I < −2B. From the above reasoning, we see that κ quantifies the contamination from the Higgsstrahlung subprocess into e + e − → ttH, and, for off-shell processes, of any additional background subprocess including contributions proportional to the HW W coupling. In table 4 we list the values of κ corresponding to the LO and NLO fits shown in fig. 9. As expected, all listed κ-values are close to 0.5. For e + e − → ttH at LO the Higgsstrahlung contribution induces a value κ > 0.5. For the off-shell process e + e − → W + W − bbH we observe a slightly larger value compared to the on-shell process, originating from additional irreducible backgrounds. The NLO QCD corrections to κ turn out to be significant. They decrease κ by 6.0% and 4.6% compared to LO for the on-and off-shell case, respectively. This can be understood from a different behavior of the signal and background contributions with respect to QCD corrections. From table 4 we can infer that at NLO interference terms are indeed negative for the on-shell ttH process. The sensitivity formula, (22), can be used to assess the impact of perturbative corrections on the extraction of y t . This is roughly half as large as the corrections reported in fig. 8. As already observed at the cross section level, the shifts in the extracted y t value that result from the inclusion of NLO corrections and off-shell contributions have comparable size and opposite sign. The magnitude of the individual effects amounts to a few per cent at 800 GeV and grows up to about 10% at full CLIC energy.

Polarization effects
We complete our study of inclusive cross sections for leptonic top-pair and Higgs associated toppair production with an investigation of possible beam polarization effects on these processes. Beam polarization is a powerful tool at linear colliders to disentangle contributing couplings and to reduce backgrounds [98] or improve the measurement of the top Yukawa coupling [13]. In tables 5 and 6 inclusive LO and NLO cross sections with different polarization settings as suggested by the favored ILC running scenarios [99] and two different collider energies are listed for the on-shell processes e + e − → tt and e + e − → ttH, respectively. While cross sections vary strongly with the beam polarization, the K-factors are unaffected. These results confirm the naive expectation that NLO QCD corrections fully factorize with respect to the beam polarization due to the uncolored initial state. On the other hand, one can view the constant K-factors in tables 5 and 6 as validation of the polarization dependent Whizard-OpenLoops-interface via a BLHA extension. The factorization also holds when top-quark decays are considered and we refrain from showing polarized cross sections for off-shell production processes.

Numerical predictions for differential distributions
Theoretically -but also experimentally -leptonic tt and ttH production and decay are very similar. Therefore, a sound understanding of tt production and decay in the continuum, where experimentally a large amount of data can easily be accumulated, is a necessary prerequisite for precision measurements of the top Yukawa coupling in the e + e − → ttH process. In this section we discuss differential predictions for e + e − → tt and e + e − → ttH at √ s = 800 GeV including NLO QCD corrections and off-shell effects in the decays. We also present predictions for the forward-backward asymmetry in e + e − → tt.

Top-pair production and decay
We start our analysis of differential distributions for top-pair production and decay considering in fig. 10 the top-quark transverse momentum distribution for the on-shell process e + e − → tt and the corresponding off-shell process e + e − → µ + ν µ e −ν e bb including leptonic decays. For the latter the top quark is reconstructed from its leptonic decay products at Monte Carlo truth level, i.e. p T,W + j b = p T, + νj b . Despite the different normalization of the two distributions, due to the fact that the on-shell process does not include leptonic branching ratios, the LO and NLO shapes are very similar below the Jacobian peak located at around 350 GeV. This peak with its large event density is smeared out by the NLO corrections, in particular due to kinematic shifts induced by the real gluon radiation, yielding corrections at the level of −20% at the peak and around +20% below the peak. For the on-shell process the phase-space above the Jacobian peak is kinematically not allowed at LO and gets only sparsely populated at NLO. In contrast, for the off-shell process this kinematic regime is allowed already at LO. The observed sizeable corrections in the transverse momentum of the intermediate top quarks also translate into relevant corrections in the directly observable transverse momentum of the final state leptons, as shown in appendix B (fig. 19). Namely, we find corrections up to −40% and up to −30% for the hardest and second hardest lepton, respectively. In a realistic setup, where experimental selection cuts have to be applied on the leptons, such effects become also relevant for the fiducial cross section in precision top physics.
Experimentally, p T,W + j b is not directly measurable, as in the considered leptonic decay mode the top quark cannot be exactly reconstructed due to the two invisibly escaping neutrinos. As a proxy we can, however, construct and measure the transverse momentum of the b-jet-lepton system, p T, + j b . Corresponding predictions for e + e − → µ + ν µ e −ν e bb are shown in fig. 11 (left). K-Factor Figure 11: Transverse momentum distribution of the bottom-jet-lepton system (left), p T, + j b , and of the j b -jb system (right), p T,bb , in e + e − → µ + ν µ e −ν e bb. Curves and bands as in Fig. 10.
Here we observe a tilt of the NLO shape with respect to the LO one, yielding corrections up to 20% for small p T, + j b and up to −40% for large p T, + j b . In contrast, the transverse momentum distribution of the j b -jb system -as shown on the right of fig. 11 -only receives mild QCD corrections at the level of 10%.
One of the observables of prime interest is the kinematic mass of the top resonance. In fig. 12 we show on the left the reconstructed invariant top-quark mass, m W + j b = m + νj b , where the + νj b system is identified based on Monte Carlo truth. At LO and close to the peak, this distribution corresponds to the Breit-Wigner that arises due to the propagator. Off-shell effects and non-resonant contributions become visible a couple of GeV away from the pole and tend to increase the background. At NLO we observe a drastic shape distortion compared to LO -in particular below the resonance peak. These NLO shape distortions are very sensitive to the cone size of the employed jet algorithm. They can be attributed to QCD radiation that escapes the b-jet forming either a separate light jet or being recombined with the other b-jet. The reconstructed invariant top-quark mass is thus on average significantly shifted compared to the top-quark resonance. Similar shape distortions have also been observed in Ref. [25] as well as at the LHC [100,101].
Again, the perfectly reconstructed top-quark mass is not directly measurable due to the escaping neutrinos. However, we can resort to the invariant mass of the b-jet and the associated charged lepton. In fact, this distribution can be used to measure the top-quark mass via [102][103][104] where m 2 j b and cos θ j b are the mean values of the corresponding invariant mass and angular distributions. Predictions for the m + j b invariant mass distribution are shown on the right of

Forward-backward asymmetries
At a future lepton collider the top quark forward-backward asymmetry A F B is defined as where θ t is the angle between the positron beam axis and the outgoing top-quark. This asymmetry can be measured with a precision below 2% [98]. The SM prediction for A F B is non-zero due to interference contributions between s-channel Z-and γ * -exchange in the dominant production process [105]. Various new physics models can substantially alter the SM prediction (for an overview cf. [106]) and thus, a precise determination of A F B serves as a stringent probe for new physics. 6 In fig. 13 we show the underlying distribution in the angle of the (reconstructed) top quark with respect to the beam axis for on-shell top-pair production and the corresponding off-shell process e + e − → µ + ν µ e −ν e bb. The prediction of a non-zero forward-backward asymmetry at lepton colliders is apparent in fig. 13 and the shape of this distribution is hardly affected by radiative corrections, which yield an almost constant K-factor of about 1.05. For cos θ W + j b 0.75, the angular distribution of the reconstructed top quark in e + e − → µ + ν µ e −ν e bb correlates with the on-shell prediction. However, for cos θ W + j b 0.75, there is an enhancement of events, which can be attributed to single-top background diagrams. This has a significant effect on the reconstructed top forward-backward asymmetry, which is reduced by about 20%, see the e + e − → µ + ν µ e −ν e bb and e + e − → W + W − bb predictions in table 7. In table 7 we list LO and NLO predictions for the forward-backward asymmetry A F B (and the corresponding asymmetry for the anti-top quark), considering different treatments of the top-quark off-shellness. In e + e − → µ + ν µ e −ν e bb, either the top-quark is reconstructed at MC truth level or the information of the neutrino momenta is dropped. NLO QCD corrections to A F B can be sizeable (up to a few percent), but are small compared to the changes associated with increasing the final-state multiplicity and taking into account all off-shell and non-resonant effects. Note that if the neutrino momenta are omitted, the relation A F B = −Ā F B is not fulfilled any more, both at LO and NLO. This can also observed directly in the angular distribution of lj b -pairs, see fig. 20 in appendix B, where there is a slightly more pronounced dip at the lower edge of cos θ l − jb than at the one of cos θ l + j b . The differences come from combinatorial issues in the event reconstruction, where the neutrino momentum in the MC truth information allows to determine the top helicity and hence the flight direction of the lepton (cf. e.g. [112]). Such information is unavailable when the neutrino kinematics are omitted.

Higgs associated top-pair production and decay
We start our analysis of differential Higgs associated top-pair production by considering in fig. 14 the energy of the Higgs boson, E H , and the invariant mass of the tt system, m tt , in the  in fig. 15. Again, we observe a strong enhancement for large Higgs boson energies and small reconstructed top-pair masses, together with a strong suppression for large reconstructed toppair masses. In contrast to the on-shell process, already at LO kinematic boundaries are washed out due to off-shell and non-resonant contributions. In particular, the E H distributions range to energies above 335 GeV, with strongly increasing NLO corrections. The m W + W − j b jb distribution at LO falls off quickly below m W + W − j b jb = 2m t , while at NLO it reaches to very small values. As already discussed in the context of fig. 12, this phase-space region is populated at NLO due to kinematic shifts of the reconstructed masses originating from the recombination of radiation from different stages of production and decay.
In fig. 16 we show the transverse momentum distribution of the reconstructed top quark and the directly observable bottom-jet-lepton system in the off-shell process e + e − → µ + ν µ e −ν e bbH. Comparing these distributions with the corresponding ones for top-pair production, shown in fig. 10 and fig. 11, we observe distinct (LO) shape differences. Instead of a pronounced peak in the p T,W + j b distribution we observe a plateau between about 100 GeV and 250 GeV. At larger transverse momenta, the distribution drops sharply to its kinematical bound at around 325 GeV. NLO QCD corrections shift both the p T,W + j b and the p T, + j b distribution towards smaller values inducing shape effects up to −50% at large p T,W + j b and up to −60% at large p T, + j b .
Finally, in fig. 17 we turn to the reconstructed kinematic top mass, m W + j b , and its directly observable relative, m + j b . We observe similar NLO shape distortions as already discussed in the case of top-pair production, shown in fig. 12. For m W + j b < m t , i.e. below the top resonance, we observe a strong NLO enhancement that translates to 20% shape corrections in the case of the m + j b distribution. As already noted before, the size of these corrections strongly depends on the details of the employed jet clustering.
Further differential distributions are shown in the appendix, cf. fig. 21.

Conclusions
In this article, we have presented the first application of the Whizard NLO framework based on a process-independent interface between Whizard and the amplitude generator OpenLoops.
We have presented a precision study for a future high-energy lepton collider considering for the first time at NLO QCD the processes e + e − → µ + ν µ e −ν e bb and e + e − → µ + ν µ e −ν e bbH, i.e. off-shell top-pair and Higgs associated top-pair production. Finite-width effects for intermediate top quarks and W bosons, single-top and non-resonant contributions as well as their interferences together with spin correlations have been taken into account consistently at NLO.
We have presented a study of inclusive cross sections varying as a function of the centerof-mass energy considering different approximations for the top off-shellness and an in-depth study at the differential level for √ s = 800 GeV. Off-shell effects play an important role even for the inclusive cross sections as the narrow-width approximation does not suffice to describe interference effects and background diagrams at energies far above threshold.
NLO QCD corrections also influence the dependence of the cross section on the top Yukawa coupling for the ttH processes, which has direct consequences for the achievable accuracy in measuring this coupling. In particular, we have shown that the NLO QCD corrections induce negative interference terms yielding a deviation from the quadratic Yukawa coupling dependence of the cross section (κ NLO < 0.5), both in the on-shell treatment of ttH production and the corresponding off-shell process.
Many facets of the strong physics potential of linear lepton colliders are based on polarized lepton beams. In order to describe beam polarization, the Binoth Les Houches Accord for the interface between Whizard and OpenLoops was generalized. As expected, we found that beam polarization has no effect on the relative size of NLO QCD corrections. It is however important to incorporate a treatment of polarization effects in the NLO framework of Whizard,  as well as QED initial state radiation (and to a lesser extent also beamstrahlung, which always factorizes), in order to allow for realistic Monte Carlo simulations in the environment of a lepton collider. In particular, as soon as one includes higher-order electroweak corrections, a consistent treatment of beam polarization is mandatory.
In addition to these inclusive studies, NLO QCD effects on differential observables have been investigated. Our results show that NLO effects can yield (clustering-dependent) corrections of up to ±50% for a variety of observables. Even larger corrections occur due to non-relativistic top-threshold effects. To obtain reliable predictions in these threshold regions, a resummed calculation is required. The easy-to-use automation of threshold matching for on-shell top quarks as well as the matching of this calculation to the relativistic continuum will be supported by Whizard in the near future. Off-shell effects are most relevant in top-mass related observables, which is of crucial importance for the determination of m t . On the other hand, Higgs observables are mostly unaffected by off-shell contributions, but can be influenced significantly by NLO QCD corrections. Studying the top-quark forward-backward asymmetry, we found the effect due to off-shell contributions dominating over NLO QCD corrections.
The study at hand was performed at the fixed-order NLO QCD level. The matching to parton showers based on an independent implementation of the POWHEG method within Whizard will be considered in the future. This will allow for realistic experimental analyses including resummation of soft and collinear radiation at the NLO+LL level. Besides the parton shower matching also NLO electroweak corrections and their interplay with the QCD corrections should be considered in the future as they are well known to play an important role for the considered processes.
In addition to its phenomenological relevance, the presented calculation demonstrates the flexibility of Whizard for NLO QCD computations at lepton colliders and the smooth interplay with OpenLoops. All distributions can be reproduced easily with these publicly available tools and finely adjusted to the experimental requirements. K-factor e + e − → µ + µ − bb Figure 18: Total cross section of the process e + e − → µ − µ + bb at LO and NLO using resonanceaware FKS subtraction. In contrast to the validation described in the text, here the physical muon mass and Higgs width have been used.
A Details of the Resonance-aware IR subtraction

A.1 Soft Mismatch
The additional soft mismatch component, which restores parts of the real-subtracted correction not covered by the resonance-aware FKS mappings, integrates the expression [34] R mism αr = dΦ B which is evaluated for each individual singular region, denoted by the index α r here. Correspondingly, R soft αr is the soft limit of the real matrix element in this region, as it is also used in the computation of real subtraction terms, and B the underlying Born matrix element. k and k res are the momenta of the radiated gluon and the intermediate resonance, respectively, and k em is the momentum of the emitter in the Born phase space. The integration is performed over the whole real phase space, which factorizes into the Born phase space dΦ B and the realradiation variables ξ, y, and φ. Note that, in contrast to the traditional FKS subtraction, where ξ = 2k 0 / √ s ≤ 1, a generalized ξ ∈ [0, ∞) is used, which originates from using integral identities. Therefore, the soft mismatch has to be evaluated with its own phase space and must be treated as a separate integration component in Whizard, additionally to Born, Real and Virtual.

A.2 Validation and Efficiency
We have checked our implementation of resonance-aware FKS subtraction extensively using the production of two massive quarks in association with two muons as a benchmark process, i.e. e + e − → bbµ + µ − . This process has only one resonance topology with two different resonance histories, Z → bb and H → bb, comprising Z pair production as well as Higgsstrahlung. We have set m b = 4.2 GeV, so that collinear divergences do not occur. For the validation of our implementation, in order to avoid any sort of cuts, we have set the muon mass equal to m µ = 20 GeV. To ensure a converging integration also in the case of the non resonance-aware approach, we have approximated the limit Γ H → ∞ by numerically fixing the Higgs width to Γ H = 1000 GeV. In this way, the standard subtraction can be compared to the improved one, see table 8, where σ real denotes the full real-subtracted matrix element and σ mism the result of the integration of the soft mismatch component. Adding the real and soft-mismatch components for the resonance-aware FKS subtraction, perfect agreement with the real radiation component of the resonance-unaware subtraction is found. Here, we want to emphasize the significantly higher number of integration calls required in the standard approach to reach the same accuracy as in the resonance-aware subtraction scheme.  Fig. 18 shows a scan of the total cross section. For this scan, we used the physical muon mass and Higgs width. There are two distinct peaks at m Z and m Z + 2m b , as well as two less pronounced enhancements at m Z + m H and 2m Z . NLO QCD corrections are in the range of +5% for √ s > 2m Z and approximately −4% for m Z + 2m b < √ s < 2m Z . Below √ s = m Z , the K-factor is significantly smaller than 1.