Scalar production and decay to top quarks including interference effects at NLO in QCD in an EFT approach

Scalar and pseudo-scalar resonances decaying to top quarks are common predictions in several scenarios beyond the standard model (SM) and are extensively searched for by LHC experiments. Challenges on the experimental side require optimising the strategy based on accurate predictions. Firstly, QCD corrections are known to be large both for the SM QCD background and for the pure signal scalar production. Secondly, leading order and approximate next-to-leading order (NLO) calculations indicate that the interference between signal and background is large and drastically changes the lineshape of the signal, from a simple peak to a peak-dip structure. Therefore, a robust prediction of this interference at NLO accuracy in QCD is necessary to ensure that higher-order corrections do not alter the lineshapes. We compute the exact NLO corrections, assuming a point-like coupling between the scalar and the gluons and consistently embedding the calculation in an effective field theory within an automated framework, and present results for a representative set of beyond the SM benchmarks. The results can be further matched to parton shower simulation, providing more realistic predictions. We find that NLO corrections are important and lead to a significant reduction of the uncertainties. We also discuss how our computation can be used to improve the predictions for physics scenarios where the gluon-scalar loop is resolved and the effective approach is less applicable.


Introduction
The top quark is the heaviest of the standard model (SM) fermions, with a mass close to the electroweak scale and a coupling to the Higgs boson close to unity.These properties could indicate its intimate connection to the dynamics of electroweak symmetry breaking (EWSB) and its role as a portal to physics beyond the SM (BSM).As an efficient topquark factory, the LHC will inspect the properties of this particle to an unprecedented level of precision, and set accurate limits on any anomalous production mechanism.The main production mechanism of top quarks at the LHC is pair production initiated by gluons, followed with roughly one third of the total cross section by single top production.
If a hypothetical new physics signal in top quark production originates from scales out of the reach of the LHC, it will manifest as low energy modifications parametrised by effective operators [1][2][3][4][5][6][7][8].Alternatively, the LHC may be able to directly produce resonances which decay into a pair of top quarks, producing a bump in a particular mass region.Resonant top-pair production has been searched for by different experiments, e.g.[9][10][11][12][13].
Resonances decaying into a t t system exist in many BSM scenarios, such as the two-Higgs-doublet model (2HDM), the supersymmetric models and several models of dynamical EWSB.Particularly relevant are the colour singlet pseudo-scalar resonances (A) H decaying dominantly to tops, pp → A/H → t t , due to the fact that in many models scalars have suppressed couplings to light fermions, proportional to their masses.In order to correctly describe this process at the LHC, and to determine the constraints on the parameters of different models, it is important to get trustworthy and robust theoretical predictions.It is now known that QCD corrections to this process have a large effect on the total cross section as well as the differential distributions, both for the SM QCD production, known at next-to-next-to-leading order in QCD and next-to-next-to-leading log soft gluon resummation [14,15] and NLO EW [16,17], and for (pseudo-)scalar production, which is known at N 3 LO in the infinite top mass limit [18]) and NLO for finite top mass [19,20].
It is also well known that the interference between resonant signal and the QCD background can drastically modify the lineshape of the signal, from a single peak to a dip-peak structure, and may be even larger than the pure signal.This effect has been studied at the leading-order (LO) in Refs.[21,22] and in more detail in Refs.[23][24][25][26].It appears to receive similar QCD corrections, as indicated by estimates in an effective field theory (EFT) approach in the soft gluon approximation [27], and in a simplified K-factor derived from the geometric mean of the background and signal K-factors [28].
In this work for the first time we provide the full perturbative calculation of the process pp(→ A/H) → t t, at NLO in QCD, including the interference, in the EFT approach, without any further approximation.The effective description is justified for BSM models in which new heavy QCD charged states provide the dominant contribution in generating gluon-scalar interactions.Our calculation does not directly apply to the cases in which the loop contribution is dominated by the top quark.However, we show that in such cases our calculation can be improved by adopting a reweighting technique, in which the Born pieces of the calculation are described by the full one-loop amplitudes [29,30].In absence of a full calculation which is currently out of reach, this approach is expected to provide a reliable estimate for the QCD corrections even in the kinematic region where the EFT approximation is not valid.Unlike Ref. [27], our reweighting includes the interference part and is done on an event-by-event basis, and can be passed to the parton shower (PS) simulation, thus providing more realistic predictions.
Our calculation is performed within the MadGraph5_aMC@NLO (MG5_aMC) [31] framework which allows us to perform automatic simulation of events with PS matching.The computation of NLO QCD corrections to top-quark processes involving anomalous cou-plings or higher-dimension operators has evolved in recent years [5,6,[32][33][34][35][36][37][38][39][40][41][42][43][44][45][46][47].In particular automated calculations within MG5_aMC have been performed in [5,6,43,46,47], and the implementation used here closely relates to these works.In the context of the EFT, one needs to consistently take into account the running and mixing between different operators, which could lead to additional complications.In the process being studied, an interesting feature that has not been discussed before is that the gluon-scalar operator mixes into the chromo-magnetic dipole-moment operator, and therefore both operators must be taken into account in a consistent NLO calculation and put together in a coherent setup.Thanks to our previous work in Ref. [5], the QCD corrections to the chromo-magnetic operator are already available.In this work, we will then deal with the remaining problems, i.e. the implementation of the operator mixing, as well as the two-loop matching which is required to determine the size of the chromo-magnetic operator.Such a calculation is an interesting one in itself: unlike the previous NLO EFT calculations in Refs.[5,6,[32][33][34][35][36][37][38][39][40][41][42][43][44][45][46][47], in this work the EFT is used in a top-down way, by starting from explicit BSM models, performing matching calculations at the same accuracy, employing RG equations down to the scale of the problem, and physical results are obtained by carrying out an NLO calculation there.
The paper is organised as follows.A general description of the features and strategy of the calculation is presented in Section 2. In Section 3 the theoretical setup is discussed in detail.In Section 4 we provide the description of several benchmark scenarios which are of phenomenological relevance.In Section 5 we show the corresponding results for each benchmark and discuss the effects of NLO corrections on the lineshape.Finally we present our conclusions in Section 6.

Heavy scalar lineshape
In this section we briefly discuss some basic features of the heavy scalar lineshape in the t t final state.In several BSM scenarios, the additional scalar couplings to fermions are hierarchical, roughly proportional to the fermion masses.This is the reason why we focus on the gluon fusion process as the dominant production mode.Contributions from q q initial states are relatively easier to deal with, and the interference between signal and background vanishes at LO for a colour singlet resonance, and so we will not consider them here.
The gluon-gluon-scalar vertex can take two forms depending on the underlying model.If the vertex is dominated by heavy particle loops (e.g.vector-like fermions or scalar top partners), or induced by strong dynamics at a high scale, this interaction will not be resolved at the scale of the scalar resonance.In this case the vertex can be simply represented by the dimension-five operators with G A µν being the gluon field strength tensor and G Aµν its dual.On the other hand, if the vertex is loop-induced with lighter particles, such as top quarks running in the loop, it will be resolved at the resonance scale and the interaction will give rise to an absorptive phase.It is convenient to distinguish between these two cases, as they often lead to different lineshapes, and the resolved case is more difficult to compute at NLO accuracy.In practice, one may need to deal with a mixed scenario, if there are contributions from both light and heavy loop particles.It has been noticed in an earlier work [22], and discussed in a series of recent works [25,26,28,[48][49][50][51][52], that the production of a heavy scalar resonance leads to large interference with the SM t t background.This large interference could be further augmented by a nontrivial relative phase between the signal and the SM background amplitude, possibly leading to more complicated structures.Possible lineshapes can vary from a pure Breit-Wigner (BW) resonance to peak-dip structures, pure dip structures, and even enhancedpeak structures, depending on the details of the underlying model [25].
To briefly explain these effects, in Figure 1 we show the loop induced resonant Feynman diagram (a), which in the heavy fermion limit can be described by a contact interaction as in diagram (b), and a SM QCD background diagram (c).Due to the large production rate in the SM, the interference is expected to be important.In particular for non-narrow resonances, the interference can be larger than the signal.
The impact of the interference on the lineshape can be understood by considering the heavy scalar propagator convoluted with the loop form factor of the top loop (we consider only the top quark loop as a resolved loop) and the contact interaction from high scale physics, which gives with ) ) where c g is the effective contact interaction coupling and c t is the normalised fermion coupling to the (pseudo-)scalar.The A H,A 1/2 loop form factors for scalar (H) and pseudoscalar (A) approach A H,A 1/2 → 1 in the limit of heavy fermion mass in the loop (τ → 0) and develop an imaginary part for a resolved loop, τ > 1.
The signal contribution comes from squaring the resonant amplitude and displays a BW shape in the m(t t) spectrum, with small perturbations from the loop form factor.The interference, on the other hand, is proportional to the real part of M sig , (2.7) In the non-resolved case, the only complex phase in this problem comes from the phase of the propagator, and thus the interference contribution depends only on the real part of the propagator, which flips sign near the resonance due to the factor of s − M 2 in its numerator.Even more interesting is the resolved case, where an additional complex phase will be provided by the loop form factors.In this case the second term in the numerator of the r.h.s of Eq. (2.7) contributes with a pure BW-shaped component to the full signal.For example, in the 2HDM, the heavy scalar and pseudoscalar couple to the gluons through top-quark loops.In this case the phase starts to appear as the mass of the resonance goes above the 2m t threshold, possibly leading to pure dip or enhanced peak structures [25].

NLO approach
Searching for heavy scalars in the t t final state is challenging, not only because of the complicated lineshapes, but also due to the top-quark invariant mass reconstruction that smears the signal, and the systematic uncertainty associated with the large production cross section of the SM top-quark pairs [51,52].For this reason, it is useful to optimise the experimental search strategy according to theoretical predictions for the lineshape, which have to include the interference effect that determines the signal shape.In fact, in a recent ATLAS search [53], the interference indeed plays an important role in the interpretation of the results in the 2HDM.On the theory side, however, the interference contribution to the signal is currently known only at the LO accuracy in QCD, while NLO corrections are expected to be important, given that the corresponding corrections to the pure signal and background are both large (at roughly ∼ 100% and ∼ 50%).Recent progress has been made in Ref. [27], where an EFT approximation together with a soft-gluon approximation has been adopted.More recently in Ref. [28], predictions based on K-factors inferred from the pure signal and background components have also been provided.
In this work, we aim to follow a more solid approach to the NLO computation of the interference, based on the EFT framework, which in the unresolved case provides accurate results (i.e.without using soft-gluon approximation), and in the resolved case, can be further improved by using reweighting techniques (by Born-reweighting or FTapprox which includes the exact real corrections [29,30]) as we will discuss in the following.These results will then be passed to PS simulation, to obtain more realistic predictions.
Let us first discuss the main difficulty at NLO.The signal is loop-induced (by light or heavy particles), therefore computing the signal at NLO requires calculating two-loop diagrams.For the signal, the initial state factorisable corrections are well known [19,20] and also implemented in Monte Carlo generators [54].Similarly the final state corrections are also well known as part of the QCD corrections to the Higgs decay width to heavy quarks [55].Non-factorisable corrections vanish for the signal, because the non-factorisable virtual corrections produce the top-quark pair in a colour octet state, and therefore do not interfere with the Born level signal diagram which contains a colour singlet.Nonfactorizable real corrections vanish for the same reason.This is however not true anymore when interference is taken into account, see Figure 2, where the interactions represented by blue dots can be either resolved or non-resolved.In this case virtual corrections that are non-factorisable could interfere with the SM QCD background, giving an NLO correction to the interference component.If the ggH(A) vertex is induced by loop particles, a complete calculation of these corrections will involve two-loop diagrams, e.g. as shown in Figure 3.These computations are difficult to perform due to the many scales involved.However, the computation of the interference in the unresolved case, mediated by a gluon Higgs effective operator, can be performed exactly at NLO, and will be presented for the first time in this work without further approximations.As we will show in the next section, including only the scalar-gluon operators (Eq.(2.1)-Eq.(2.2)) is not enough for a consistent calculation in the EFT.This is because operator mixing into a top-quark dipole moment becomes relevant in this process, and we will have to include this operator in our approach.This procedure also involves the matching stage, where a two-loop matching needs to be performed also for the dipole operator.
We note here that the exact results in the EFT can be also used to improve the predictions for the resolved case.It is well known that the EFT approach does not hold if s goes above 4m 2 X , where X is the particle running in the loop, but reweighting with the exact Born-level result helps to improve the predictions for heavy scalar single or double production.In particular it has been shown that the EFT approximation with a Born reweighting is a good approximation for single Higgs production even for heavy scalar mass [56] and thus we follow a similar approach to improve our predictions for the resolved case.While this is appropriate for the signal, for the interference a simple reweighting based on the ratio between LO exact and EFT results is not reliable, since reweighting is rather problematic as the ratio I exact /I EFT diverges as the exact and EFT amplitudes cross zero at different points.In order to address this problem a phase reproducing the exact LO behaviour can be introduced to the EFT amplitude, by adding an ad-hoc imaginary part to the effective operator coefficient.

Theoretical setup
In this section we discuss the technical setup of our approach in more detail.

EFT for the unresolved case
We first discuss the case where the ggH(A) vertex is mediated by high-scale physics.In this case the EFT is a good approximation.Consider the following EFT of the SM augmented by a heavy scalar H where O HG is given in Eq. (2.1).
If the full theory is known, one can determine the coefficient by explicitly computing the vertex in the full theory.For example, if this contribution is induced by a vector-like fermion, F , via the following Yukawa interaction: by computing diagrams shown in Figure 4 one finds C HG (µ) By using this operator one can compute the process at NLO in QCD, which involves at most one-loop diagrams.For consistency, corrections to Eq. (3.3) from two-loop matching should be incorporated [19]: However, this computation leads to uncancelled UV poles from non-factorisable contributions, whenever a top quark, a gluon and a scalar form a loop, as shown in Figure 5.
The reason for the UV divergence is not difficult to understand.The effective Lagrangian, given in Eq. (3.1), is in principle not a renormalizable one.It is well known that the EFT is renormalisable only if one considers the complete set of higher-dimensional operators up to a certain dimension.However, since so far only the operator O HG is added, it could mix into a different operator not included in the Lagrangian, leading to non-renormalisability of the theory.In fact, through the triangles shown in Figure 5, the O HG induces a chromo-magnetic dipole operator of the top quark, The physical interpretation is the following.The right-hand side of the second diagram in Figure 5 is a triangle loop made by a Higgs, a top, and a gluon.In the full theory, it represents a similar diagram but with the blue dot replaced by a fermion loop, namely the Barr-Zee diagram [57], an important contribution for lepton dipole moments from additional scalars.We see here that the same diagram is playing a role also for the top quark, leading to a chromo-magnetic dipole moment of the top.Two comments are in order.First, the gg → h calculation in the SM is based on essentially the same Lagrangian, but the renormalisability is not an issue there, because O tG simply does not enter the gg → h process.Here we see that the problem occurs because we have a scalar that decays into the t t final state, which makes O tG relevant.This also explains why the problem occurs only in the non-factorisable pieces.Second, one option to avoid this problem is to simply assume that the gluons in these loops are always soft for the dominant contribution, which is the assumption used in Ref. [27].The assumption is a good one for the resonant region, but its validity for predicting the lineshape in a much larger range is not guaranteed.The process being considered here is most interesting if the interference effect is large, which implies a large width and a non-trivial lineshape over a much larger range of the m(t t) distribution.For this reason, in this work we will not use this approximation.
It is clear that the solution is to take Eq.(3.1) and extend it to incorporate also the O tG operator and perform the calculation.With the correct mixing coefficient, proper counterterms will be generated to cancel the UV poles of the diagrams in Figure 5.
The missing piece now is to determine the coefficient of C tG in the EFT via a matching procedure.Without the matching, the contribution from Figure 5 will depend on the scheme in which we determine the counterterms.

Matching
At the one loop level there is no contribution to C tG .However to be consistent with the NLO calculation, we need to consider the matching at two loops, just like what we have done for C HG in Eq. (3.4).This can be done by computing the Barr-Zee contribution to the gt t function, and matching the result to the EFT contribution.The actual computation can be model dependent.For a model with a CP-even scalar H coupled to vector-like fermions F through a Yukawa coupling given by Eq. (3.2), this is illustrated in Figure 6.On the left-hand side, the full result for the Barr-Zee diagrams has been computed in the context of lepton dipole moments.We take the analytical expression from Ref. [58] and expand in powers of 1/m F .On the right-hand side, the one-loop contribution from O HG is computed using C HG which has been matched at one loop.The counterterm from C tG is determined in the M S scheme.Setting Λ = m F , and neglecting the top-quark mass, the matching is represented by the following equation: The left-hand side is finite as expected from the full theory, and involves two scales, m F and m H .The same log m H dependence appears also on the right-hand side, with exactly the same coefficient.This is expected because the two theories are supposed to describe the same IR physics.On the other hand, the loop contribution from the right-hand side is divergent, also expected because the EFT modifies the UV structure of the theory.By choosing µ = m F , we can cancel the log terms in the matching, and the resulting coefficients are defined at the scale m F .The UV pole and the remaining finite terms from the loop contribution on the right have only tree-level structure.The former will be cancelled by the last term, i.e. an M S counterterm, while the latter will match the full result on the left-handed side by choosing a proper value for C tG .We find: This together with Eq. (3.4) defines our EFT at the two-loop accuracy.Scheme dependence always cancel in the final results as long as the matching and the actual calculation based on EFT are performed within the same scheme.This is because the sum of the last two terms in Eq. (3.9) is scheme independent.The top-quark mass, m t , has been neglected in the matching.This will only give rise to power corrections such as (m t /m H ) 2 .Given the smallness of C tG , this contribution is expected to be negligible.
Adding a CP-odd scalar, A, will not affect the above approach.Consider extending Eq. (3.2) by where A is a CP-odd scalar.The theory is still CP-conserving, and after matching will lead to our complete EFT: with operators written explicitly in Eq. (2.1), Eq. (2.2) and Eq.(3.5).Note that C tG is real as long as CP is conserved.The matching can be performed in a similar way since the two-loop corrections for both ggA and ttg (from a CP-odd scalar) are known [19,59].We find It is not our purpose to present matching results for all BSM extensions with additional scalars.However, the procedure outlined above can be carried out for any perturbative BSM model.

Running
The actual calculation and simulation will be performed at the renormalisation scale µ, conventionally chosen at µ = m H /2 (m A /2).In case the scale m F and µ are well separated, one needs to run the EFT from m F down to µ to resum the large logarithmic contributions.
For O(α s ) mixing this is done by solving the RG equations where for the three operator coefficients C HG , C A G, C tG , the matrix γ ij is given by (3.17) The coefficients at a given scale µ are thus given by where β 0 = 11 − 2/3n f , and n f = 5 is the number of running flavors.
The above procedures take into account the O(α s ) mixing from O HG,A G to O tG .The mixing from O tG to O HG also exists, but is of order O(y 2 t ), and is negligible because We also note here that unlike the matching, the running in EFT is model-independent.

Calculation and automation
Once C i (µ) is known, NLO predictions based on the EFT can be computed.In this work we use the automated framework based on MG5_aMC.Calculations at NLO with higherdimensional operators have been recently performed for the top-quark sector of the SMEFT [5,6,43,46,47].For example, in Ref. [47] predictions for pp → t th are obtained at NLO in QCD for O tG and the SM Higgs O HG operators.The approach has the advantage that results are fully exclusive, automatically matched to PS through MC@NLO [60], and can be directly used in experimental analyses.
A similar implementation including O tG and O HG has been created for this study.There are a few differences compared with previous works: • The SM Higgs field is replaced by heavy scalar(s). 1 The O HG (and O A G) is defined with a prefactor of g 2 s , so that the LO signal and background have the same power of α s .The relevant mixing is thus from O HG to O tG , not the other way around like in pp → t th.
• Complex mass scheme [61,62] is used in order to take into account the widths of heavy scalars.
The framework is, however, very similar and has been fully tested in previous works.

Resolved case and reweighting
The approach described above provides the exact NLO QCD prediction for the EFT, which is a good approximation if the ggH (ggA) vertex is induced by heavy particles with masses larger than m H(A) /2.For lighter particles, in particular the top quark which always couples directly to the scalar, the EFT is not valid anymore.In addition, it fails to capture the absorptive part of the loop, which is crucial for understanding the lineshape.Because of the phase, a simple reweighting does not provide a solution.In practice, the EFT and top-loop amplitudes become zero at different phase-space points, therefore the amplitude ratio |M exact | 2 /|M EFT | 2 used for Born-reweighting can diverge.As already discussed in the previous section, a possible but ad hoc solution consists of introducing a phase in the EFT amplitude, mimicking the absorptive part of the loop.This is achieved by adding an imaginary part to the operator coefficient in the following simple form The values of a and b are constant and are obtained at Born level by a fitting procedure that ensures that the exact and EFT amplitudes cross zero at the same mass and with the same gradient.In more detail, the resolved amplitude is computed using the exact top-mass dependence in the form factor A H,A 1/2 of Eq. 2.3.The corresponding EFT result obtained for the infinite top mass limit, i.e. setting A H,A 1/2 to 1, is multiplied at the amplitude level by (a + bi).These two forms of the amplitude are then used to compute the interference contribution to the partonic cross-section for gg → t t.This partonic cross section is then examined as a function of the invariant mass of the top-antitop pair.For the resolved amplitude, we numerically extract the invariant mass at which the interference is zero, and also compute the slope of the partonic cross-section at the same point.Requesting that the EFT result crosses zero at the same invariant mass, and with the same slope provides a system of two equations with two unknowns which we solve to extract the value of a and b.This procedure ensures the ratio |M exact | 2 /|M EF T | 2 remains finite at all phase space points and allows us to perform an event-by-event reweighting [63].The reweighting leads to the exact result at LO and a Born-improved result at NLO.We note here that with our setup we can also include the exact real emission amplitudes in a fashion similar to what has been done for double Higgs production in [29,30], but we refrain from doing so for simplicity.A discussion of the real emission amplitudes in the interference for this process can be found in [28].

Benchmarks
To study the lineshape, we will perform the calculation for several benchmark models, which cover both the resolved and unresolved cases.The impact of QCD corrections in our EFT approach will be illustrated by the first two benchmark models, where the ggH(A) vertex is not resolved.In the first scenario, we consider a CP-even scalar that couples to a heavy vector-like fermion doublet, which induces the ggH vertex at one loop.The matching procedure discussed in the previous section can be explicitly carried out.In the second model, we consider a CP-odd scalar as a pseudo-Goldstone boson of new strong dynamics [64,65], whose decay into top-quarks may shed light on an eventual dynamical fermions mass generation [66,67].It is important to mention that in these two benchmarks the top-loop contributions also exist and as they are resolved cannot be described by the EFT.For this reason these two benchmarks should be interpreted as ideal cases that help us to understand the impact of NLO corrections in a scenario where the EFT approximation holds.To make them more realistic, we require that the top-loop contributions are always subdominant.
For our final benchmark we consider the CP-even and CP-odd Higgses of the 2HDM.In this case the top-loop induced ggH vertex is fully resolved as only top and bottom quarks run in the loops.For these benchmarks we construct the EFT and improve our predictions by Born-reweighting.
In the following we give more details for the three benchmark models.

Benchmark A
Consider a vector-like quark doublet, F , with Yukawa coupling to a CP-even scalar H described by Eq. (3.2).We choose the following parameters: The ggH vertex from the F fermion running in the loop is given by Eq. 3.4, multiplied by a factor of two.For the top-quark loop, we simply replace m F by m t and y F by y t / √ 2. This does not fully capture the top-loop contribution, in particular the phase, however given that the dominant contribution is coming from the F fermion, the EFT is still a very good approximation.To demonstrate this, we plot the partonic cross sections computed from the EFT and from the exact one-loop amplitude in where we have taken m t to be 172.5 GeV.The above values define our EFT.

Benchmark B
As a benchmark of CP-odd state we consider the models of partial compositeness [68] of Ref. [67], more specifically the a state of M3 model with (n ψ , n χ ) = (−4, 2) in the η decoupling limit, α = ζ.We choose the compositeness scale f = 800 GeV and the mass m a = 1 TeV.The relevant couplings of A ≡ a in Eq. (3.12) are 2 : The O A G operator is generated through the axial-vector anomaly which couples the technipions to a pair of gauge bosons.We assume The total width is dominated by the top decay and has some correction from di-gluon decay and is given by Γ A = 37.5 GeV .(4.8) It is interesting to note that C A G and ỹt have the same sign and this implies that the top loop has an opposite contribution to the gluon coupling, creating a dip-peak structure instead of the usual peak and dip.This opposite contribution also generates some strong cancellations if one approximates the full top loop with a heavy fermion, as we have done for benchmark A, thus for this benchmark we use only the effective coefficient generated by the high scale physics and neglect the top loop in the gluon coupling.This is illustrated in Figure 8, where the partonic cross section normalised to the QCD background computed at LO is displayed.The difference between the curves stems from the way to treat the resolved top loop in the amplitude in Eq. (2.3).The exact curve has the full form factor structure.The heavy fermion limit m t → ∞ and A A 1/2 (τ ) → 1 is a bad approximation because of the large cancellation between c g and c t .A better approximation, given by the blue curve, is to use the limit where m t → 0 and the top loop vanishes, which is justified by the heavy scales we are probing. 2The conversion from parameters defined in Ref. [67] to our convention in Eq. (3.12) is given by: Type-I and II Table 1.Top quark Yukawa couplings to the light (heavy) CP-even and CP-odd Higgs bosons.These are identical for type-I and type-II 2HDM.

Benchmarks C1 and C2
For our final BSM scenario we employ the 2HDM [69], which introduces a second SU (2) L doublet Φ 2 and gives rise to five physical Higgs bosons: one light (heavy) neutral, CP-even state h (H); one neutral, CP-odd state A; and two charged Higgs bosons H ± .In this work we take h to be the 125 GeV Higgs.The input parameters determining all properties of a 2HDM scenario are: with the convention 0 ≤ β − α < π (with 0 < β < π/2).
The ratio of the Yukawa couplings over the SM values are shown in Table 1 for both type-I and type-II 2HDMs.Electroweak precision tests, the LHC Higgs results and searches for heavy scalar particles, along with unitarity, perturbativity and vacuum stability constrain the parameter space of the model.In the selection of 2HDM benchmarks, these constraints are taken into account included through the public tools 2HDMC [70], Hig-gsBounds [71,72] and HiggsSignals [73,74].The two benchmarks we will employ are shown in Table 2 Where we have taken m t to be 172.5 GeV.The above values define our EFT.For this scenario, the EFT is not expected to accurately describe the lineshape as the top-loop is resolved.The ratio over the SM background is shown at LO in Figure 9 for the interference (I) and signal (S).In particular the interference lineshape differs between the one-loop and EFT calculation, due to the phase difference between the amplitudes.A solution is to introduce a complex phase, in the ad-hoc form of C AG Λ → (a + bi) × C AG Λ , where a and b are extracted by matching the exact and EFT amplitudes at leading order in particular at the point where they cross zero.For this benchmark introducing this phase can be used to improve the description of the lineshape as shown in Figure 9.This is particularly useful also for the NLO computation which is improved by Born reweighting.

C2: CP-even and CP-odd resonances
The corresponding couplings and widths for this benchmark are:  The above values define our EFT.Similarly to benchmark C1 the EFT does not provide a reliable prediction of the lineshape and a phase (one for the scalar and one for the pseudoscalar resonance) can be used in order to match the lineshape of the interference.The results are shown in Figure 10.

Results
In this section we present results for the four benchmarks described in section 4. We fix our EFT scale at one half of the scalar resonance mass, while for renormalisation and factorisation scales we choose a dynamical scale equal to one half of the sum of the transverse masses in the final state.We vary the scales independently by a factor of two up and down to estimate the scale uncertainties.The Yukawa coupling between the H and the top quark is renormalised by the M S scheme but fixed at the scale m H /2. We use the LO and NLO NNPDF2.3 parton distribution functions [75] for the corresponding computation.We obtain results both at fixed-order and also matched to the parton shower with MC@NLO.
For the parton shower we always use Pythia8 [76].The fully inclusive phase-space is considered and the tops are kept stable, even though their decays can be considered in a fully automatic way via the MadSpin [77] package allowing further analysis of exclusive and realistic observables.The cross sections for the SM background, the interference, and the signal components at the LHC 13 TeV are given in Table 3.Large QCD corrections can be observed for the signal, but not for the interference.This however does not mean that the radiative corrections are not important for the interference, as the interference displays a cancellation below and above the resonance.In fact, looking at the scale uncertainties, we can see that our calculation improves the precision by a factor of ∼ 4.

Benchmark
To further investigate our results, in Figure 11 we show the fixed order m t t distribution.Signal (S) and interference (I) lineshapes are displayed separately.The signal peaks at the resonant mass 500 GeV, and the radiative correction is quite large as expected.In particular, we observe an increasing K-factor as the t t invariant mass moves away from the resonance to lower values.This is mainly because the heavy Higgs decays to t t + g, where the gluon reduces the total energy of the t t pair, so that events near the peak at LO may be shifted to the left at NLO.A similar but less significant effect can be observed also on the interference lineshape, which shows a typical peak-dip structure.In the lower panels, we find that radiative corrections to the rate are large, while the improvements on the scale uncertainties are mild.Finally, as a comparison, we follow the approach in Ref. [28] and rescale the LO lineshape by a K-factor inferred by using √ K S K B , where K S and K B are the (bin by bin) K-factors for the SM and the signal.As shown in the plot, the resulting lineshape is a reasonable approximation in the absence of a full NLO computation for the interference which we provide in this work, but it tends to overestimate the size of the interference at and below the resonance.In addition, by construction this approach does not improve the LO scale uncertainty, as can be observed in the last panel.
In Figure 12, the same results are shown but matched to the PS.Compared with Figure 11, the size of the radiative corrections is now tamed below the resonance, mainly because the emission from decay products is partly taken into account by the PS.The Kfactor for the signal is in general flat, with a small enhancement near the peak.In contrast, the K-factor for the interference shows a peak-dip structure, mainly because the zero point has been slightly shifted.This effect is not reproduced by a K-factor rescaling used in Ref. [28].As one can see in the last panel, the inferred K-factor √ K S K B is flat and does not reproduce the peak-dip structure, and as a result this approximation will predict an interference lineshape that is shifted to the left compared with our full calculation.Finally, compared with the fixed-order results, the showered lineshape has a slightly lower peak due to smearing effects.
In Figure 13, we show the total BSM effect, including signal and interference.In particular, we compare the sum of NLO signal and LO interference with the full NLO prediction, and their ratio is given in the lower panel.The comparison shows the impact of this work.The rate below the resonance is increased by about 50%, while the increase above the resonance is even larger.A shift of the zero point is also observed.Again, we compare with the K-factor rescaling approximation used in Ref. [28].Similar to Figure 12, we find that this approximation predicts a lower lineshape, i.e. it tends to underestimate the peak, but overestimate the dip.It also underestimates the resonant mass of the scalar.In addition, in the same plot we show the sum of LO signal and LO interference, normalised to SM LO background.The full NLO results normalised to NLO background gives a much higher peak, which is expected because the QCD correction to the resonance is larger than that to the background.
In the same graph we also plot the error due to the PS treatment of the interference.The origin of the error is explained as follows.In general, the PS algorithm as implemented in Pythia8 distinguishes between events with and without an explicit resonant state.If the resonance is present, showering is treated in a factorised way.External coloured legs of the production part are generated, in such a way so that the invariant mass of the resonance is conserved.If an intermediate resonance is not present, then the PS only sees the full process, and external coloured legs will be generated without conserving the invariant mass of the resonance.Another difference, particular to this process, is that events with and without a colour neutral resonance correspond to different classical colour flows.Both facts affect the results of showering.In practice, we can treat the SM background as non-resonant while the signal part as resonant.However, whether the interference events should be considered with a resonant intermediate state is not well-defined, and this is the origin of what we call PS uncertainty of the interference.The problem is not so severe for processes with small interference contributions.However, the process we consider here is a special one, where it is the interference contribution that really determines the shape of the BSM contribution.Therefore it is important to assess the possible consequences of this uncertainty.
By default, samples generated by MG5_aMC contain both resonant and non-resonant events, depending on the relative size of signal and background.Such a scheme is not a physical one, as it implies that the shape of the interference will be determined by the size of the effective operator coefficient.To estimate the uncertainty due to this fact, we generate two additional sets of events, one with only resonant events, the other with only non-resonant.We then pass these two sets of events through the PS, and the resulting m t t distributions form an error band, which is displayed in Figure 13 by a light blue band.We can see that this effect is not important away from the resonance, however it does shift the lineshape near the resonant area, where the interference is expected to be large.The normalised uncertainty is given in the lower panel by the red band.Again we see that in general this uncertainty is small, except for the near-resonance region, where this effect can become comparable to the scale uncertainties.
In Refs [78] and [79,80] resonance aware matching has been proposed within the MC@NLO and Powheg frameworks respectively, with focus on narrow coloured resonances.In our configuration the resonances are broader and the interference plays a more important role.In the case of the large interference the ambiguities related to the treatment of the resonance are more severe, and our estimate of the parton shower uncertainty could well be an underestimate.
In general the PS uncertainty arises from the intrinsic limitations of the current PS algorithms that are publicly implemented, and therefore a solution to this problem is beyond the scope of this work.We estimate its potential impact to benchmark A, but we do not expect a significant difference for the other benchmarks.For this reason we will not consider this error for the rest of the paper.The cross section for Benchmark B is shown in Table (4) with a generation cut applied to the invariant mass of the top pair system, m(t t) > 750 GeV.This is to ensure that we have enough statistics in the high m(t t) region near the resonance.

Benchmark
In Figure 14 we show the invariant mass distribution of the top pair system m t t at fixed order in perturbation theory.In the upper panel the NLO and LO predictions are shown normalised bin-by-bin to the background QCD prediction.Pure signal (S) and interference (I) are shown separately.Results obtained using the approximate average K-factor, √ K S K B , are also shown.The lower panels show the K-factors and the scale uncertainty as an envelope.For the signal, the same behaviour observed for the scalar benchmark can be noticed: a large effect at the lower tail of the distribution can be observed as the gluon emission reduces the energy of the t t system.In the high energy tail a smaller scale uncertainty is obtained with the NLO calculation.For the interference, the QCD corrections increase the rates mainly below the resonance for the fixed-order predictions.It is interesting to restate that in this scenario there is a dip-peak structure due to the opposite sign contribution of the high energy strong sector.The approximate K-factor, √ K S K B , overestimates the size of the corrections, in particular in the region below the resonance.
In Figure 15 we show the equivalent distribution with matching to the PS.The conclusions are similar to those described in the previous benchmark.at the resonance mass, while the interference K-factor is mostly flat for the NLO+PS predictions.The overall K-factors are however lower than in benchmark A, due to the higher mass considered.The cross sections for the signal, background and their interference for scenarios C1 and C2 are shown in Table 5   the central value of the prediction but also reduce the scale uncertainties.Differential results for the 2HDM benchmarks C1 and C2 are shown in Figures 16  and 17 respectively.The invariant mass distribution of the top quark pair is computed at (N)LO+PS, with the signal (S) and interference (I) contributions and their corresponding K-factors displayed separately.The LO results are exact, i.e. they include the full top mass dependence, while for the NLO results we use the phase-improved EFT to compute the NLO corrections and then apply a Born-reweighting at an event-by-event level, as discussed in the previous section.As the reweighting is based on event generation we present only (N)LO+PS for these benchmarks.Results obtained using the approximate average Kfactor, √ K S K B , are also shown.For both benchmarks the signal K-factors peak around the resonant masses.The interference K−factor shows a peak-dip structure as the NLO corrections shift the zero crossing point from its LO position.The approximate K-factor, √ K S K B is much flatter, with small peaks driven by the peaks in the signal K-factor at the resonance masses.Moreover, the uncertainties of that prediction are much larger as they are by construction LO.

Benchmarks C1 and C2
For benchmark C1 we also show in yellow the prediction obtained for the interference in the EFT limit, i.e. without any reweighting.This result should be close to the approximation in [27].As expected from the partonic results shown in Figure 9 the EFT prediction does not a provide a good description of the interference lineshape.
We stress that for scenarios where two resonances are present the experimental resolution will be crucial to establish whether one or two peaks are visible in the invariant mass distribution.For our benchmark the mass difference of 150 GeV should be sufficiently large.We note that here the interference distribution suffers from low statistics in the region between the two resonance masses, where the cross-section is very small.

Conclusions
We have computed for the first time the NLO QCD corrections for the interference between signal and background in resonant scalar or pseudoscalar top pair production using an EFT approach.The interference between signal and background is crucial in the determination of the resonant lineshape in a series of motivated BSM scenarios, where new heavy scalar and/or pseudoscalar particles are coupled to the top quarks.In the case where the gluon fusion scalar production is dominated by heavy particles running in the loop, or strong dynamics introducing a point-like interaction between the scalar and the gluons, the EFT provides a very good approximation of the process, which we have verified by analytical calculation.The computation at NLO in the EFT limit reduces the problem to one-loop level and can be performed without further approximations.
Using this approach, we have computed the NLO corrections within the MG5_aMC framework, to obtain an accurate description of the t t lineshape.As a nontrivial feature of this calculation, we found that the effective gluon-scalar operator mixes into the chromomagnetic dipole operator, which implies that both operators, together with their mixing effects, need to be incorporated in the computation to guarantee a finite and physical result.To demonstrate how the coefficients of both operators can be extracted from the underlying theory, we have performed two-loop matching to a UV complete model, and evolved the operator coefficients to the scale where the calculation will be performed.Whilst the matching is model dependent, the running and mixing between the operators is model independent.Furthermore, to handle the cases where the EFT is less applicable, i.e. when there are light particles running in the gluon fusion loop, we have demonstrated that the EFT results can still be used to improve the predictions.This can be done by introducing a complex phase in the coefficient of our effective operators to match the absorptive part of the one-loop amplitude, and employing a Born-reweighting using the exact LO amplitudes to obtain results beyond LO.
Our setup is fully automated and allows us to compute results both at the fixed order and matched to the parton shower.We have presented results for a representative set of benchmarks, covering both scalar and pseudoscalar resonances in the unresolved and resolved cases.In particular, we have studied a scenario with a scalar resonance coupling to a heavy vector-like quark doublet, a scenario of a pseudoscalar state in a model of partial compositeness, and two benchmarks in the 2HDM.The first two scenarios can be well described by the EFT, whilst for the 2HDM scenarios, a complex phase and reweighting have been used to improve the EFT predictions.In all cases we have found that QCD corrections are important and significantly reduce the scale uncertainties of the predictions.At the differential level the corrections lead to nontrivial K-factors and therefore are crucial for a precise determination of the lineshape.We examined both fixed-order and NLO+PS distributions, and in particular assessed the impact of the uncertainty due to the treatment of the interference contribution in the matching to the PS simulation.We found that this uncertainty is relevant in the region close to the resonance where it can become comparable to the scale uncertainties.Finally, we compared our predictions with previous approximate NLO results, which are based on the geometric mean of the respective signal and background K-factors.We found significant improvements in the resonant region.
In summary, we have computed the NLO QCD corrections to resonant scalar production and decay into top quarks in the EFT limit taking into account the signal and background interference.Our computation improves the accuracy and precision of the theory predictions, allowing for more reliable analyses of resonant top pair production.Matching to the PS simulations is provided in an automated way, and thus predictions can be used in realistic simulations in the context of resonant searches at the LHC.Combined with advances in the experimental techniques, we expect that our work will improve the sensitivity of experimental searches for BSM scalars, by reducing the systematical uncertainty from the theory side, and by optimising the experimental strategies according to an accurate lineshape description.On the theory side, it will also allow us to extract reliable constraints on the parameters of various BSM models.Finally, the calculation itself is an interesting one, as the two-loop matching to both operators allows for the EFT to be used in a top-down way, to improve the theoretical predictions in a specific BSM scenario.

2. 1 Figure 1 .
Figure 1.LO signal and background.Only one diagram from SM background is shown.

Figure 3 .
Figure 3. Two-loop contributions involved in non-factorisable corrections.Loop particles depend on the details of the model.

Figure 4 .
Figure 4. One-loop matching to compute C HG .

Figure 6 .
Figure 6.Two-loop matching to determine the coefficient C tG at two-loop accuracy.

πFigure 8 .
Figure 8. Partonic cross section for gg → t t, including the resonance signal and the interference, normalised to the t t background, computed from EFT and from the exact one-loop amplitude.

Figure 9 . 11 × 10 − 8
Figure 9. Signal (S) and interference (I) over background ratios, computed from EFT (with and without the additional phase) and from the exact one-loop amplitude.

Figure 10 .
Figure 10.Signal and interference over background ratios, computed from EFT (with and without the additional phase) and from the exact one-loop amplitude.

Figure 11 .
Figure 11.Signal and interference lineshapes at fixed order (N)LO for LHC 13 TeV.Results are normalised to the SM fixed order at NLO (fNLO) lineshape.Lower panels show the K-factors and the scale uncertainties of signal and interference respectively.

Figure 12 .
Figure 12.Signal and interference lineshapes at (N)LO matched with PS for LHC 13 TeV.Results are normalised to the SM fNLO lineshape.Lower panels show the K-factors and the scale uncertainties of signal and interference respectively.

Figure 13 .
Figure 13.Total BSM effects (signal+interference) with the interference at LO/NLO, at the LHC 13 TeV matched to PS. Results are normalised to the SM fNLO lineshape, except for the LO curve which is normalised to SM fLO.The lower panel displays the K-factors and the scale uncertainties of signal and interference respectively.

Figure 14 .
Figure 14.Signal and interference lineshapes at fixed order (N)LO for LHC 13 TeV.Results are normalised to the SM fNLO lineshape.Lower panels show the K-factors and the scale uncertainties of signal and interference respectively.

Figure 16 .
Figure 16.Signal and interference lineshapes at (N)LO matched with PS for LHC 13 TeV.Results are normalised to the SM fNLO lineshape.Lower panels show the K-factors and the scale uncertainties of signal and interference respectively.

Figure 17 .
Figure 17.Signal and interference lineshapes at (N)LO matched with PS for LHC 13 TeV.Results are normalised to the SM fNLO lineshape.Lower panels show the K-factors and the scale uncertainties of signal and interference respectively.
Partonic cross section for gg → t t at LO, including the resonance signal and the interference, normalised to the t t background, computed from EFT and from the exact one-loop amplitude.

Table 2 .
and are briefly discussed below.Type tan β sin(β − α) m H m A m H ± Parameter choices for the 2HDM benchmarks used in our study.All masses are given in GeV.The light Higgs mass is fixed to m h = 125 GeV.The heavy scalar width is negligible while the pseudoscalar width is Γ A = 7.35 GeV and the top quark decay branching fraction is approximately 65%.Following Eqs.(3.4), (3.10), and (3.18), we find the following operator coefficients, at scale µ EF T = m A /2:

Table 3 .
Cross sections for benchmark A in pb.Uncertainties are from renormalisation and factorisation scale variation.

Table 4 .
Cross sections for benchmark B in pb with m(t t) > 750 GeV.Uncertainties are from renormalisation and factorisation scale variation.
The signal K-factors peaks

Table 5 .
Cross sections for benchmark C1 in pb.Uncertainties are from renormalisation and factorisation scale variation.
and 6 respectively.The QCD corrections change significantly Signal and interference lineshapes at (N)LO matched with PS for LHC 13 TeV.Results are normalised to the SM fNLO lineshape.Lower panels show the K-factors and the scale uncertainties of signal and interference respectively.

Table 6 .
Cross sections for benchmark C2 in pb.Uncertainties are from renormalisation and factorisation scale variation.