Electroweak top-quark pair production at the LHC with $Z'$ bosons to NLO QCD in POWHEG

We present the calculation of the NLO QCD corrections to the electroweak production of top-antitop pairs at the CERN LHC in the presence of a new neutral gauge boson. The corrections are implemented in the parton shower Monte Carlo program POWHEG. Standard Model (SM) and new physics interference effects are properly taken into account. QED singularities, first appearing at this order, are consistently subtracted. Numerical results are presented for SM and $Z'$ total cross sections and distributions in invariant mass, transverse momentum, azimuthal angle and rapidity of the top-quark pair. The remaining theoretical uncertainty from scale and PDF variations is estimated, and the potential of the charge asymmetry to distinguish between new physics models is investigated for the Sequential SM and a leptophobic topcolor model.


Introduction
The Standard Model (SM) of particle physics is a very successful theory describing a wealth of experimental data up to collision energies of 13 TeV reached at CERN's Large Hadron Collider (LHC). This includes the recent observation of a Higgs-like particle with a mass of 125 GeV that seems to corroborate the simplest description of electroweak symmetry breaking [1][2][3]. However, the SM is based on the unintuitive semi-simple gauge group SU(3) C ×SU(2) L ×U(1) Y , that together with the running behavior of the associated gauge couplings intriguingly points towards a larger unification at some higher mass scale. The simple gauge group SU (5) can accomodate the complete SM gauge group and its 15 fermions, but not a right-handed neutrino, and it is in addition strongly disfavored by searches for proton decay. It also does not allow to restore parity symmetry and does not provide a natural solution to the neutrino mass hierarchy. Both of these important and perhaps related problems are solved in simple gauge groups of higher rank like E 6 or SO (10), that can be broken consecutively as in E 6 →SO(10)×U(1) ψ and SO(10)→SU(5)×U(1) χ , respectively. Parity restoration is achieved in left-right symmetric models, SU(3) C ×SU(2) L ×SU(2) R ×U(1) Y , which together with other models of similar group structure, but different quantum number assignments form a class of general lowerscale models, commonly called G(221) models. They have recently been classified [4], and their phenomenology has been studied not only at the LHC [5][6][7], but also in ultrahighenergy cosmic rays [8]. Common to all these possible extensions of the SM is their prediction of a new heavy neutral gauge boson (Z ), that is associated with the additional SU(2) or U(1) subgroup after symmetry breaking [9,10]. In many cases, the Z boson can decay leptonically, making it a prime object of experimental searches at the LHC. For simplification, these searches are mostly based on the (theoretically unmotivated) Sequential SM (SSM), where the Z boson couples to other SM particles like the SM Z boson. In this model and the leptonic (i.e. Drell-Yan) channel, the ATLAS and CMS collaborations have already excluded Z bosons with masses below 2.90 TeV [11] and 2.96 TeV [12], respectively. For a recent overview of experimental mass limits see Ref. [6], where it is also shown that for certain G(221) models the mass limits are enhanced to 3.2-4.0 TeV, when higher-order QCD corrections are included.
In this paper, we focus not only on the SSM, but also on a situation where the Z boson does not couple to leptons, but preferentially to top quarks, so that the above mass limits are invalidated. Models of the G(221) class, where processes of the Drell-Yan type are inaccessible at the LHC, include leptophobic (LP), hadrophobic (HP) and fermiophobic (FP) models, whereas left-right (LR), un-unified (UU) and non-universal (NU) models remain accessible. The LP model with a W -boson mass of about 2 TeV has been put forward as a possible explanation for the excesses of W Z and W h production observed recently by ATLAS and CMS at the LHC [13]. As the heaviest particle in the SM with a mass of 173 GeV [14], the top quark may very well play a special role in electroweak symmetry breaking. This motivates, e.g., the NU model, where the first and second SU(2) gauge groups couple exclusively to the first/second and third generation fermions, respectively. It also motivates models with new strong dynamics such as the topcolor model [15,16], which can generate a large top-quark mass through the formation of a top-quark condensate. This is achieved by introducing a second strong SU(3) gauge group which couples preferentially to the third generation, while the original SU(3) gauge group couples only to the first and second generations. To block the formation of a bottom-quark condensate, a new U(1) gauge group and associated Z boson are introduced. Different couplings of the Z boson to the three fermion generations then define different variants of the model [17]. A popular choice with the LHC collaborations is the leptophobic topcolor model (also called Model IV in the reference cited above) [18], where the Z couples only to the first and third generations of quarks and has no significant couplings to leptons, but an experimentally accessible cross section.
The strongest limits on Z bosons arise of course from their Drell-Yan like decays into electrons and muons at the LHC. This is due to the easily identifiable experimental signatures [6]. The top-pair signature is more difficult, as top quarks decay to W bosons and bottom quarks, where the latter must be tagged and the two W bosons may decay hadronically, i.e. to jets, or leptonically, i.e. into electrons or muons and missing energy carried away by a neutrino. In addition and in contrast to the Drell-Yan process, the electroweak top-pair production cross section obtains QCD corrections not only in the initial, but also in the final state. For conclusive analyses, precision calculations are therefore extremely important to reduce theoretical uncertainties, arising from variations of the renormalization and factorization scales µ r and µ f and of the parton density functions (PDFs) f a/p (x a , µ f ), and for an accurate description of the possible experimental signal and the SM backgrounds.
At the LHC, the hadronic top-pair production cross section obtains up to next-to-leading order (NLO) the contributions where the numerical indices represent the powers of the strong coupling α S (µ r ) and of the electromagnetic coupling α, respectively. The first and third terms representing the SM QCD background processes qq, gg → tt and their NLO QCD corrections, including the qg channel, have been computed in the late 1980 [19][20][21][22]. Furthermore, NLO predictions for heavy quark correlations have been presented in [23], and the spin correlations between the top quark and antiquark have been studied in the early 2000s [24,25]. The fourth term represents the electroweak corrections to the QCD backgrounds, for which a gauge-invariant subset was first investigated neglecting the interferences between QCD and electroweak interactions arising from box-diagram topologies and pure photonic contributions [26] and later including also additional Higgs boson contributions arising in 2-Higgs doublet models (2HDMs) [27]. The rest of the electroweak corrections was calculated in a subsequent series of papers and included also Z-gluon interference effects and QED corrections with real and virtual photons [28][29][30][31][32]. In this paper, we focus on the second and fifth terms in Eq. (1.2) (highlighted in red), i.e. the contribution σ 0;2 for the Z signal and its interferences with the photon and SM Z boson and the corresponding QCD corrections σ 1;2 . Due to the resonance of the Z boson, we expect these terms to be the most relevant for new physics searches. A particular advantage of this choice is that the calculation of σ 1;2 can then be carried out in a model-independent way as long as the Z couplings are kept general, whereas the fourth term σ 2;1 is highly model-dependent due to the rich structure of the scalar sector in many models. The sixth term in Eq. (1.2) is suppressed by a relative factor α/α s with respect to the fifth and thus small. The production of Z bosons (and Kaluza-Klein gravitons) decaying to top pairs has been computed previously in NLO QCD by Gao et al. in a factorized approach, i.e. neglecting all SM interferences and quark-gluon initiated diagrams with the Z boson in the t-channel, and for purely vector-and/or axial-vector-like couplings as those of the SSM [33]. We have verified that we can reproduce their K-factors (i.e. the ratio of NLO over LO predictions) of 1.2 to 1.4 (depending on the Z mass) up to 2%, if we reduce our calculation to their theoretical set-up and employ their input parameters. Their result has triggered the Tevatron and LHC collaborations to routinely use a K-factor of 1.3 in their experimental analyses (see below). The factorized calculation by Gao et al. has been confirmed previously in an independent NLO QCD calculation by Caola et al. [34]. Like us, these last authors include also the additional quark-gluon initiated processes and show that after kinematic cuts they reduce the K-factor by about 5 %. However, they still do not include the additional SM interferences, which they claim to be small for large Z -boson masses. As we will show, this is not always true due to logarithmically enhanced QED contributions from initial photons. In contrast to us, they also include top-quark decays in the narrow-width approximation with spin correlations and box-diagram corrections to interferences of the electroweak and QCD Born processes (σ 2;1 in Eq. (1.2)), which are, however, only relevant for very broad resonances. If the (factorizable) QCD corrections to the top-quark decay are included, the K-factor is reduced by an additional 15%. The globally smaller K-factor of Caola et al. is thus explained by calculational aspects and not by different choices of input parameters.
The SM backgrounds are today routinely calculated not just in NLO QCD, but at NLO combined with parton showers (PS), e.g. within the framework of MC@NLO or POWHEG [35,36]. A particularly useful tool is the POWHEG BOX, in which new processes can be implemented once the spin-and color-correlated Born amplitudes along with their virtual and real NLO QCD corrections are known and where the regions of singular radiation are then automatically determined [37]. Calculations of this type have already been performed by us in the past for the Drell-Yan like production of Z bosons [38], heavy-quark production in the ALICE experiment [39], and the associated production of top quarks and charged Higgs bosons [40,41]. In this work, we provide a calculation of the Z signal with a final top-quark pair at the same level of accuracy, including all interferences with SM Z bosons and photons as well as the logarithmically enhanced QED contributions from initial-state photons, which we will discuss in some detail. We also present details about the spin-and color-correlated Born amplitudes, the treatment of γ 5 and renormalization procedure in our calculation of the virtual corrections, as well as the validation of our NLO+PS calculation, which we have performed with the calculation for Z bosons of Gao et al. at NLO [33] and for tree-level and one-loop SM matrix elements with MadGraph5 aMC@NLO [42] and GoSam [43].
Experimental searches for resonant top-antitop production have been performed at the Tevatron and at the LHC mostly for the leptophobic topcolor model with a Z -boson coupling only to first and third generation quarks [17,18]. In this model, the LO cross section is controlled by three parameters: the ratio of the two U(1) coupling constants, cot θ H , which should be large to enhance the condensation of top quarks, but not bottom quarks, and which also controls both the Z production cross section and decay width, as well as the relative strengths f 1 and f 2 of the couplings of right-handed up-and down-type quarks with respect to those of the left-handed quarks. The LO cross sections for this model are usually computed for a fixed small Z width, Γ Z = 1.2% × m Z , effectively setting the parameter cot θ H , and the choices f 1 = 1, f 2 = 0, which maximize the fraction of Z bosons that decay into top-quark pairs. We have verified that we can reproduce the LO numerical results in the paper by Harris and Jain [18] for Z masses above 1 TeV and relative widths of 1% and 1.2%, but not 10%, if we neglect all SM interferences. As stated above, the LO cross sections are routinely multiplied by the experimental collaborations by a K-factor of 1.3 [13]. At the Tevatron with center-of-mass energy √ S = 1.96 TeV and in the lepton+jets top-quark decay channel, CDF and D0 exclude Z bosons with masses up to 0.915 TeV [44] and 0.835 TeV [45], respectively. The weaker D0 limit can be explained by the fact that CDF use the full integrated luminosity of 9.45 fb −1 , while D0 analyze only 5.3 fb −1 and furthermore do not use a K-factor for the signal cross section. At the LHC, the ATLAS and CMS collaborations have analyzed 20.3 fb −1 and 19.7 fb −1 of integrated luminosity of the √ S = 8 TeV LHC run employing the K-factor of 1.3. The result is that narrow leptophobic topcolor Z bosons are excluded below masses of 1.8 TeV and 2.4 TeV, respectively [46,47]. At the LHC, the CMS limit is currently considerably stronger than the one by ATLAS despite the slightly smaller exploited luminosity. The reason is that CMS performed a combined analysis of all top-quark decay channels (dilepton, lepton+jets and all hadronic), while ATLAS analyzed only the lepton+jets channel. For Γ Z = 10% × m Z , the CMS mass limit is even stronger and is found to be 2.9 TeV. We emphasize that the narrow width assumption employed in most experimental analyses need not be realized in nature and that in this case a proper treatment of SM interference terms as provided in our full calculation is required.
The LHC has just resumed running with an increased center-of-mass energy of 13 TeV, which is planned to be increased to 14 TeV in the near future. We therefore provide numerical predictions in this paper for both of these energies and for two benchmark models, i.e. the SSM and the leptophobic topcolor model. The predictions for the SSM are readily obtained by taking over the Z -boson couplings from the SM, with the consequence of again a relatively small width Γ Z 3% × m Z for Z masses between 3 and 6 TeV. We focus on the invariant-mass distribution of the top-quark pair, which is the main observable exploited for resonance (and in particular Z -boson) searches, but also show results for the distributions that are most sensitive to soft parton radiation beyond NLO, i.e. the transverse momentum p tt of the top-antitop pair and their relative azimuthal angle φ tt .
The forward-backward asymmetry A F B of top-antitop events with positive vs. negative rapidity difference between the two has also been suggested as a very useful observable to distinguish among different models [48]. At the Tevatron (a pp collider, where top quarks are produced predominantly in the direction of the proton beam), long-standing discrepancies of CDF and D0 measurements with the SM prediction at NLO [49,50] have triggered numerous suggestions of new physics contributions [48], e.g. of light Z bosons coupling in a flavor non-diagonal way to up and top quarks [51]. Only recently the SM prediction at next-to-next-to-leading order (NNLO) [52] has been brought in agreement with the newest inclusive measurement by CDF [53] and differential measurement by D0 [54]. At the LHC (a pp collider), a charge asymmetry A C can be defined with respect to the difference in absolute value of the top and antitop rapidities [55]. We therefore also provide numerical predictions for this observable in our two benchmark models and at current and future LHC center-of-mass energies.
Our paper is organized as follows: In Sec. 2 we present analytical results of our calculations at LO and the NLO virtual and real corrections, including details about SM interference terms, our treatment of γ 5 , our renormalization procedure and the subtraction method employed for the soft and collinear divergences in the real corrections. In Sec. 3 we discuss the implementation of our calculation in POWHEG and present in particular the color-and spin-correlated Born amplitudes, the definition of the finite remainder of the virtual corrections, the implementation of the real corrections with a focus on the rather involved treatment of QED divergences, and the validation of our tree-level matrix elements in the SM against those of the automated tool MadGraph5 aMC@NLO [42] and of the virtual corrections against those of GoSam [43] as well as of our numerical pure Z -boson results against those obtained by Gao et al. and Caola et al. Our new numerical predictions for the LHC are shown and discussed in Sec. 4, and Sec. 5 contains our conclusions. Several technical details of our calculation can be found in the Appendix.

NLO QCD corrections to electroweak top-pair production
In this section, we present in detail our calculation of the NLO QCD corrections to electroweak top-pair production through photons, SM Z bosons and additional Z bosons with generic vector and axial-vector couplings to the SM fermions. We generate all Feynman diagrams automatically with QGRAF [56] and translate them into amplitudes using DIANA [57]. The traces of the summed and squared amplitudes with all interferences are then calculated in the Feynman gauge and D = 4 − 2ε dimensions in order to regularize the ultraviolet (UV) and infrared (IR) divergences using FORM [58]. Traces involving the Dirac matrix γ 5 are treated in the Larin prescription [59] by replacing γ µ γ 5 = i 1 3! ε µνρσ γ ν γ ρ γ σ . To restore the Ward identities and thus preserve gauge invariance at one loop, we perform an additional finite renormalization for vertices involving γ 5 .

Leading-order contributions
The leading-order (LO) Feynman diagrams contributing to the electroweak production of top-quark pairs at O(α) through photons, SM Z bosons and new Z bosons are shown sum-marily in Fig. 1. The cross section dσ/dt, differential in the Mandelstam variable t denoting the squared momentum transfer, is then obtained by summing all three corresponding amplitudes, squaring them, summing/averaging them over final-/initial-state spins and colors and multiplying them with the flux factor 1/(2s) of the incoming and the differential phase space 1/(8πs) of the outgoing particles. The Mandelstam variable s denotes the squared partonic center-of-mass energy. The result, given here for brevity only in four and not D dimensions, is where B qq is the modulus squared of the Born amplitude averaged/summed over initial/final spins and colors, V, V ∈ {γ, Z, Z }, the superscript q denotes the flavor of the incoming massless quarks, s, t, u are the partonic Mandelstam variables, and m t is the top-quark mass. Note that we use the Pauli metric, in which the dot-product has an overall minus sign with respect to the Bjorken-Drell metric [60]. The terms D V , D V stem from the propagator denominators and take the usual form To take into account the finite widths of the Z and Z bosons, we have introduced complex are proportional to the axial (A) and vector (B) couplings of the various gauge bosons to the massless quarks (q = u, d, s, c, b) and the top quark (t), where s W (c W ) are the sine (cosine) of the weak mixing angle θ W , Q q is the fractional charge of quark flavor q, and a q V and b q V are the model-dependent vector and axial-vector couplings of the Z and Z bosons, e.g. a u for all up-and down-type quarks in the SM. Although individual interference terms may contain imaginary parts, they cancel as expected after summation.

One-loop virtual corrections
The one-loop virtual corrections contributing to electroweak top-pair production at O(α s α 2 ) originate from the interferences among the one-loop diagrams shown in Fig. 2 with the treelevel diagrams in Fig. 1. Note that one-loop electroweak corrections to the QCD process qq → g * → tt have zero interference with the electroweak diagrams in Fig. 1, since such contributions are proportional to the vanishing color trace Tr(T a ). In particular, the interference term of the box diagram in Fig. 3 with the amplitudes in Fig. 1  it would of course contribute at O(α 2 s α). As already mentioned, the virtual amplitudes are regularized dimensionally. The appearing 30 distinct loop integrals are then reduced to a basis of three master integrals using integration-by parts identities [61,62] in the form of the Laporta algorithm [63] as implemented in the public tool REDUZE [64,65]. One is thus left with the evaluation of three master integrals: the massive tadpole, the equal-masses two-point function, and the massless two-point function. The solutions of these integrals are well known [66]. For completeness, we provide their analytic expressions in App. A.

vanishes, whereas
In dimensional regularization, the UV and IR singularities in the virtual corrections appear as poles of 1/ε and 1/ε 2 . Since neither the couplings nor the top-quark mass have to be renormalized at NLO, the UV singularities can be removed by simply adding the Born cross section multiplied with the quark wave-function renormalization constants ψ∈{q,q,t,t} We use the on-shell renormalization scheme, in which δZ q = 0 for the initial-state massless quarks and for the final-state top quarks. Since we are using the Larin prescription for γ 5 (see above), we must perform an additional finite renormalization to restore the Ward identities. The corresponding constant has been calculated up to three loops in the MS scheme [59]. At one loop, it reads and multiplies all appearing factors of γ 5 . Once the UV divergences are renormalized, we are left with infrared collinear and soft divergences that match the correct structure given for instance in Refs. [67,68]. For completeness, we provide the analytic expressions of the IR poles in App. B.

Real emission corrections
At O(α S α 2 ), the following 2 → 3 tree-level processes contribute: (i) q +q → t +t + g and (ii) g + q(q) → t +t + q(q).  As a consequence of the KLN theorem, the soft and soft-collinear divergences cancel in the sum of the real and virtual cross sections, while the collinear singularities are absorbed into the parton distribution functions (PDFs) by means of the mass factorization procedure. The singularities in the real corrections are removed in the numerical phase space integration by subtracting the corresponding unintegrated counter terms [67,68]. The fact that the collinear divergences appearing in Figs. 5 (c) and (d) involve a photon propagator has two consequences: (i) we have to introduce a PDF for the photon inside the proton and (ii) the corresponding underlying Born process shown in Fig. 6, g + γ → t +t, must be included in the calculation. The squared modulus of the corresponding Born amplitude, Figure 6. Photon-induced top-pair production of O(α S α). These diagrams must be added for a consistent subtraction of the collinear singularities.
averaged/summed over initial/final state spins and colors, is with Q t the fractional electric charge of the top quark (2/3), Although this process is formally of O(α S α) and thus contributes to σ 1;1 , it is multiplied by a photon distribution inside the proton of O(α), so that the hadronic subprocess p + p → g + γ → t +t is effectively of O(α S α 2 ). As we will see in Sec. 4, this channel is indeed numerically important.

POWHEG implementation
We now turn to the implementation of our NLO corrections to electroweak top-pair production, described in the previous section, in the NLO+PS program POWHEG [37]. We thus combine the NLO precision of our analytical calculation with the flexibility of parton shower Monte Carlo programs like PYTHIA [69] or HERWIG [70] that are indispensible tools to describe complex multi-parton final states, their hadronization, and particle decays at the LHC. Since the leading emission is generated both at NLO and with the PS, the overlap must be subtracted, which is achieved using the POWHEG method [36] implemented in the POWHEG BOX [37]. In the following, we describe the required colorand spin-correlated Born amplitudes, the definition and implementation of the finite remainder of the virtual corrections, and the real corrections with a focus on the subtleties associated with the encountered QED divergences. All other aspects such as lists of the flavor structure of the Born and real-emission processes, the Born phase space, and the four-dimensional real-emission squared matrix elements have either already been discussed above or are trivial to obtain following the POWHEG instructions [37]. We end this section with a description of the numerical validation of our implementation.

Color-correlated Born amplitudes
The automated calculation of the subtraction terms in POWHEG requires the knowledge of the color correlations between all pairs of external legs i, j. The color-correlated squared Born amplitude B ij is formally defined by for two incoming (i, j = q,q) or outgoing (i, j = t,t) particles and zero otherwise.
As we have seen in Sec. 2.3, we also have to include the gluon-photon induced pair production process in order to treat the QED divergence occurring in the gq real-emission correction. We thus also have to calculate the color-correlated squared Born matrix element for this process. The color structure of the corresponding Feynman diagrams, see Fig. 6, factorizes in the amplitude, and we can thus directly calculate the color-correlated in terms of the averaged/summed modulus squared of the Born matrix element with color factor C = N C C F = (N 2 C −1)/2. Applying Eq. (3.1) to all pairs of colored external legs, we obtain As is easily verified, a completeness relation coming from color conservation holds: and similarly for B 41 + B 43 . These cross checks are also performed automatically in POWHEG.

Spin-correlated Born amplitudes
The spin-correlated squared Born amplitude B µν j only differs from zero, if leg j is a gluon. It is obtained by leaving uncontracted the polarization indices of this leg, i.e.
where M({i}, s j ) is the Born amplitude, {i} represents collectively all remaining spins and colors of the incoming and outgoing particles, and s j is the spin of particle j. The polarization vectors ε µ s j are normalized according to Similarly to the color-correlated Born amplitudes, we have a closure relation, namely where B is the squared Born amplitude after summing over all polarizations. Since processes without external gluons lead to vanishing contributions, we must only consider the gluon-photon induced top-pair production and then modify POWHEG in such a way that the subtraction terms for the QED divergence in the gq channel can also be constructed. We therefore compute here explicitly the expression for B µν 2 , where the subscript 2 designates the photon leg (see Fig. 6). Applying the above procedure then leads to where 14) As for the color-correlated squared Born matrix element, the closure relation of Eq. (3.9) is implemented in POWHEG as a consistency check.

Implementation of the virtual corrections
For the implementation in POWHEG, the virtual corrections must be put into the form (3.19) General expressions for the coefficients a and c ij can be found, e.g., in App. B of Ref. [71] and in Refs. [72,73]. µ r is the renormalization scale, and Q is an arbitrary scale first introduced by Ellis and Sexton [74] and identified in POWHEG with µ r . The finite part V fin. is then obtained form our calculation of the virtual corrections in Sec. 2.2.

Real corrections and QED divergences
Like the Born contributions, the real-emission squared amplitudes have been implemented in POWHEG for each individual flavor structure contributing to the real cross section. As already stated above, the diagram in Fig. 5 (a) is finite and does not involve any singular regions. The diagrams in Fig. 4 and Fig. 5 (b) have the same underlying Born structure as the LO process qq → tt, followed or preceded by singular QCD splittings of quarks into quarks (and gluons) or of gluons into quarks (and antiquarks), so that their singular regions are automatically identified by POWHEG. The diagrams in Fig. 5 (c) and (d) involve, however, the photon-induced underlying Born diagrams in Fig. 6, preceded by a singular QED splitting of a quark into a photon (and a quark). The corresponding QED singularities were so far not treated properly in POWHEG. Only the singular emission of final-state photons had previously been implemented in Version 2 of the POWHEG BOX in the context of the production of single W bosons [75] and the neutral-current Drell-Yan process [76].
We therefore also implemented the photon-induced Born structures in Fig. 6, replaced the POWHEG subtraction for the QCD splitting of initial quarks into gluons (and quarks), which doesn't occur in our calculation, by a similar procedure for the QED splitting of initial quarks into photons (and quarks), and enabled in addition the POWHEG flag for real photon emission, which then allows for the automatic factorization of the initial-state QED singularity and the use of photonic parton densities in the proton. Note that this also restricts the possible choices of PDF parametrizations, as photon PDFs are provided in very few global fits.

Validation
Our implementation of the electroweak top-pair production with new gauge-boson contributions has been added to the list of POWHEG processes under the name PBZp. It allows for maximal flexibility with respect to the choices of included interferences between SM photons and Z bosons as well as Z bosons, the vector and axial-vector couplings of the latter, and the choices of renormalization and factorization scales (fixed or running with p 2 T + m 2 t or s) in addition to the standard POWHEG options. The SM Born, real and 1/ε-expansion of the virtual matrix elements have been checked against those provided by MadGraph5 aMC@NLO [42] and GoSam [43], respectively. After including the Z -boson contributions, we checked our full implementation with respect to the cancellation of UV and IR divergences. We validated, in addition to the renormalization procedure described in Sec. 2.2, the completeness relations for the color-and spin-correlated Born amplitudes and performed the automated POWHEG checks of the kinematic limits of the real-emission amplitudes. In particular, we have checked explicitly that the variable describing the collinear QED singularity shows a regular behavior after the implementation of our new QED subtraction procedure. Restricting ourselves again to the SM, our total hadronic cross section with the qq initial state only could be shown to fully agree with the results in MadGraph5 aMC@NLO, which does not allow for a proper treatment of the QED divergence in the gq initial state.
As already discussed in the introduction, the production of Z bosons decaying to top pairs has been computed previously in NLO QCD by Gao et al. in a factorized approach for purely vector-and/or axial-vector-like couplings as those of the SSM [33]. They neglected, however, all SM interferences and quark-gluon initiated diagrams with the Z boson in the t-channel. We can reproduce their K-factors of 1.2 to 1.4 (depending on the Z mass) up to 2%, if we reduce our calculation to their theoretical set-up and employ their input parameters. In the independent NLO QCD calculation by Caola et al. [34], the authors include also the additional quark-gluon initiated processes and show that they reduce the K-factor by about 5 %. However, they still do not include the additional SM interferences, which they claim to be small for large Z -boson masses. As we have discussed in detail, this is not always true due to the logarithmically enhanced QED contributions from initial photons. If we exclude SM interferences and the (factorizable) QCD corrections to the top-quark decay, we can also reproduce their K-factors.

Numerical results
In this section, we present numerical results for electroweak top-quark pair production including Z -boson contributions at LO and NLO from our new POWHEG code [37], which we coupled to the parton shower and hadronization procedure in PYTHIA 8 [69]. Our results pertain to pp collisions at the LHC with its current center-of-mass energy of √ S = 13 TeV. Only for total cross sections, we also study how much the reach in Z mass is extended in a future run at √ S = 14 TeV. The top quark is assigned a mass of m t = 172.5 GeV as in the most recent ATLAS searches for Z bosons in this channel [46] and is assumed to be properly reconstructed from its decay products. At the top-pair production threshold, α(2m t ) = 1/126.89. The values of sin 2 θ W = 0.23116, m Z = 91.1876 GeV and Γ Z = 2.4952 GeV were taken from the Particle Data Group [10]. The width of the Z boson depends on its mass and its sequential Standard Model (SSM) or leptophobic topcolor (TC) couplings. We vary the mass for total cross sections between 2 and 6 TeV and fix it to 3 TeV for differential distributions. As stated in Sec. 1, in the case of TC the Z width is set to 1.2% of its mass, and the couplings are f 1 = 1 and f 2 = 0. We use the NNPDF23 nlo as 0118 qed set of parton densities fitted with α s (m Z ) = 0.118, which includes the required photon PDF and allows to estimate the PDF uncertainty [77,78]. The renormalization and factorization scales are varied by individual factors of two, but excluding relative factors of four, around the central value µ r = µ f = √ s. In contrast to the two existing NLO calculations [33,34], which take only the Z -boson exchange and no SM interferences into account and where m Z was chosen as the central scale, our choice of √ s also applies to the SM channels and interpolates between the different physical scales appearing in the process.

Total cross sections
To illustrate the total number of events to be expected from resonant-only Z -boson production at the LHC, we show in Fig. 7  of √ S = 13 TeV in the SSM (dashed red curve) and TC (dashed black curve), together with the associated renormalization and factorization scale uncertainties (blue bands) and PDF uncertainties (green bands). As one can see, in the case of the SSM (lower curves) the PDF uncertainty is larger than the scale uncertainty in the entire range of m Z masses from 2 to 6 TeV considered here. Conversely, for the TC model (upper curves), it is the scale uncertainty which dominates for m Z 5 TeV, while the PDF uncertainty takes over only at larger values of m Z , since the PDFs at large momentum fractions x a,b are less precisely known. The uncertainties at NLO (note that the PS don't affect the total cross sections) are about ±15% at low masses and increase to ±35% in the SSM and ±20% in TC at higher masses. For an integrated luminosity of 100 fb −1 , the number of expected events falls from 10 4 for m Z = 2 TeV to 10 for m Z = 6 TeV in the SSM and is about   an order of magnitude larger in TC. When the LHC energy is increased to 14 TeV, the corresponding total cross sections (full curves) at high Z -boson mass are larger by about 50%, and the mass reach is extended by about 500 GeV, less of course than the increase in the hadronic energy √ S, of which only a fraction is transferred to the initial partons and the hard scattering.
Even for resonant-only Z -boson production, the K-factor is not completely massindependent, as can be seen in Fig. 8. In TC (lower plot), it increases only modestly from 1.3 to 1.45 in the mass range considered here, while in the SSM (upper plot) it increases much more from about 1.45 to 1.85. In contrast, it depends very little on the LHC centerof-mass energy of 13 TeV (open circles) or 14 TeV (full circles). In this figure, the scale and PDF uncertainties can also be read off more precisely than in the previous figure.
In Tab. 1 we list the total cross sections in LO for top-pair production at O(α 2 s ), O(α s α) and O(α 2 ) in the SM, SSM and TC, i.e. including the SM backgrounds, together with the corresponding NLO corrections. The Z -boson mass is set here to 3 TeV, and for our LO predictions we use the NNPDF23 lo as 0119 qed PDF set, since a set with α s (m Z ) = 0.118 is not available at this order. Comparing first the LO results only, we observe that the pure QCD processes of O(α 2 s ) have a total cross section of about 474 pb, i.e. two orders  (1) of magnitude larger than the photon-gluon induced processes of O(α s α) with 4.87 pb as naively expected from the ratio of strong and electromagnetic coupling constants in the hard scattering and in the PDFs. The suppression of the pure electroweak with respect to QCD processes is more than three orders of magnitude, as expected from the ratio of coupling constants in the hard scattering and when taking into account that the QCD processes have both quark-and gluon-initiated contributions. The Z -mediated processes in the SSM and TC have only cross sections of 5 and 12 fb, respectively compared to 366 fb from the SM channels alone, which therefore clearly dominate the total electroweak cross sections. The interference effects are destructive in the SSM (−4%), but constructive in TC (+2%).
When a cut on the invariant mass of the top-quark pair of 3/4 of the Z mass (i.e. at 2.25 TeV) is introduced, the SM backgrounds are reduced by more than three orders of magnitude, while the signal cross sections drop only by about 10%. The interference effects then become more important in the SSM (−7%), but not in TC (+2%) with its very narrow Z width of 1.2% of its mass. While an invariant-mass cut strongly enhances the signal-to-background ratio, the LHC experiments still have to cope with signals that reach only 3 to 8 % of the QCD background, which makes additional cuts on kinetic variables necessary.
The NLO corrections for the QCD processes are well-known and can be computed with the published version of POWHEG (HVQ) [79]. At the LHC with its high gluon luminosity, the qg channels opening up at NLO are known to introduce large K-factors, here of about a factor of three. The NLO corrections for the purely electroweak processes are new even in the SM, where we have introduced a proper subtraction procedure for the photon-induced processes. The K-factors for the qq channel are moderate in the SM (+56%), SSM (+58%) and TC (+56%), where the last two numbers are dominated by SM contributions and therefore very similar. Only after the invariant-mass cut the differences in the models become more apparent in the K-factors for the SM (±0%), SSM (+19%) and TC (+23%). However, similarly to the QCD case the qg channel, and also the γg channel opening up for the first time at this order, introduce contributions much larger than the underlying Drell-Yan type Born process. Note that the LO γg cross section computed with NLO α s and PDFs must still be added to the full NLO qq +gg cross sections. An invariantmass cut is then very instrumental to bring down the K-factors and enhance perturbative stability, as one can see from the LO γg and in particular the NLO results in the SSM and TC.

Differential distributions
We now turn to differential cross sections for the electroweak production of top-quark pairs that includes the contribution of a SSM or TC Z boson with a fixed mass of 3 TeV.
The invariant-mass distributions of top-quark pairs in Fig. 9 exhibit steeply falling spectra from the SM background from 10 −2 to 10 −7 pb/GeV together with clearly visible resonance peaks of SSM (top) and TC (bottom) Z bosons at 3 TeV, whose heights and widths differ of course due to the different couplings to SM particles in these two models. In particular, the TC resonance cross section is about an order of magnitude larger than the one in the SSM in accordance with the total cross section results in the previous subsection (see Fig. 7). What becomes also clear from the lower panels in Fig. 9 (top and bottom) is that the K-factors are highly dependent on the invariant-mass region and can reach large factors around the resonance region. This is particularly true for TC (bottom), but also for the SSM, and related to the fact that the position of the resonance peak is shifted towards lower invariant masses from LO to NLO due to additional radiation at this order. As one can see, this effect is already present if parton showers are added to the LO calculation, so that the NLO+PS to LO+PS comparison mostly results in an increased K-factor at and above the resonance.
The effect of interferences between SM and new physics contributions is shown in Fig.  10, where the sum of the squared individual contributions (blue) is compared with the square of the sum of all contributions (green) in the SSM (top) and TC (bottom). As one can see, the interference effects shift the resonance peaks to smaller masses, and their sizes are reduced. When the ratios of the two predictions are taken (lower panels), it becomes clear that predictions without interferences overestimate the true signal by a factor of two or more.
The two variables that are particularly sensitive to soft-parton radiation and the associated resummation in NLO+PS Monte Carlo programs are the net transverse momentum of the observed particle (here top-quark) pair (p tt ) and the azimuthal opening angle between them (φ tt ), which are 0 and π, respectively, at LO. At NLO they are balanced by just one additional parton and thus diverge and exhibit physical behavior and turnover only at NLO+PS, i.e. after resummation of the left-over kinematical singularities. These well-known facts can also be observed in Figs. 11 and 12, where for obvious reasons the LO δ-distributions at 0 and π are not shown. As expected, the NLO (green) predictions diverge close to these end points, while the NLO+PS (red) predictions approach finite asymptotic values. Again, a similar behavior is already observed at LO+PS accuracy, although with different normalization and shape. Interestingly, the resummation works much better for purely Z -mediated processes (lower panels) than if SM and interference contributions are included (upper panels). This effect can be traced back to the fact that in the SM-dominated full cross section the top-pair production threshold at 2m t = 345 GeV is almost one order of magnitude smaller than the mass m Z = 3 TeV governing the exclusive Z -boson channel.
In our discussion of total cross sections in Sec. 4.1, we had included analyses of scale and PDF uncertainties at NLO, but not of the uncertainty coming from different PS implementations, as the PS does not influence total cross sections, but only differential distributions. To estimate this uncertainty, we therefore show in Figs. 9 and 11 also results obtained with the HERWIG 6 PS (dashed red) [70] in addition to those obtained with our standard PYTHIA 8 PS (full red) [69]. The dashed red curves in the lower panels of Fig.  9 represent the ratios of the HERWIG 6 over the PYTHIA 8 PS results. As one can see there, the invariant-mass distributions in the SSM and TC are enhanced by the HERWIG 6 PS at the resonance at 3 TeV by about 10%, while the region just below it is depleted by a smaller amount, but over a larger mass region. The PS differences are therefore smaller (by factors of three to six, except for the PDF error in TC) than those of the scale and PDF uncertainties in Fig. 8. The SSM transverse-momentum distribution in Fig. 11 falls off a bit faster with the HERWIG 6 PS than with the PYTHIA 8 PS at large transverse momenta, while in TC it is slightly enhanced at low values, but no significant differences appear between the angularly ordered HERWIG 6 PS and the dipole PS in PYTHIA 8.
The importance of next-to-leading-logarithmic (NLL) contributions that go beyond the leading-logarithmic (LL) PS accuracy can be estimated by a comparison with analytic NLL resummation calculations. These have not been performed for top-quark, but only for lepton final states [38]. In Fig. 5 of this paper, it has been found that the invariant-mass distribution shows no significant difference, while the LL transverse-momentum distribution computed with the HERWIG 6 PS is somewhat smaller than the one obtained with NLL resummation, but that it stays within the residual scale uncertainty of the latter.
Rapidity distributions of the top-quark pair are shown in Figs. 13 and 14. If SM contributions are taken into account (top), they are much flatter than if only the heavy resonance contributes (bottom), i.e. the top-quark pairs are then produced much more centrally. The effect is similar, but somewhat less pronounced in TC (Fig. 14) than in the SSM (Fig. 13) due to the broader resonance in this model. Even for rapidity distributions NLO effects are not simply parametrizable by a global K-factor, as it varies from 1.6 to 2.1, when SM contributions are taken into account (blue curves in the upper K-factor panels) and drops from 1.6 to 1.4 or even below, if they are not taken into account (blue curves in the lower K-factor panels). As expected, the parton showers (green curves in the K-factor panels) have little effect on the central parts of the rapidity distributions, and they only slightly influence the forward/backward regions through additional parton radiation from the initial state.
A particularly sensitive observable for the distinction of new physics models is the forward-backward asymmetry defined at pp colliders, where ∆y = y t − yt is the rapidity difference of top and antitop quarks, and the somewhat more complex charge asymmetry defined at pp colliders, where ∆|y| = |y t | − |yt| is the corresponding difference in absolute rapidity [55]. In Fig. 15, the sensitivity of A C to distinguish between the SSM (top) and TC (bottom) is confirmed, as this observable exhibits very different magnitudes at the resonance (11 ± 1% vs. ±0.1%) and far below it (2.5 ± 0.5% in both plots), where the SM contributions dominate. Since A C is defined as a ratio of cross sections, NLO and PS corrections cancel to a large extent and are barely visible above the statistical noise. Only for TC, where the rapidity distribution in Fig. 14 (lowest panel) showed distinct features in the ratio of NLO+PS/LO+PS, the transition from the low-mass to the resonance region happens more abruptly in fixed order (NLO) than with PS. If we assume an integrated luminosity of 100 fb −1 and integrate over an invariant-mass window of 100 GeV around the resonance peak at 3 TeV, one would expect 10 −5 pb/GeV×100 fb −1 × 100 GeV = 100 events. A 10% asymmetry in the SSM then implies a difference of 10 events with an error of 3, so that A C = (10 ± 3)%. This would be sufficient to distinguish the SSM from the SM and TC.

Conclusions
In this paper we presented the calculation of the O(α S α 2 ) corrections to the electroweak production of top-antitop pairs through SM photons, Z and Z bosons, as predicted in the Sequential SM or in tecnicolor models. Our corrections are implemented in the NLO parton shower Monte Carlo program POWHEG. Z reconances are actively searched for by the ATLAS and CMS experiments at the LHC with its now increased center-of-mass energy of 13 TeV. We have consistently included interferences between SM and new physics contributions and have introduced a proper subtraction formalism for QED singularities. With a great variety of numerical predictions, we have demonstrated the mass dependence of the K-factor, the changing relative sizes of scale and PDF uncertainties, the large impact of new partonic channels opening up at NLO (in particular of those induced by photon PDFs in the proton), and the non-negligibility of interference effects. Distributions in invariant mass were shown to be particularly sensitive to the latter. The all-order resummation of perturbative corrections implicit in the parton shower has been shown to make the transverse-momentum and azimuthal angle distributions of the top-antitop pair finite and physical. Heavy new gauge-boson contributions were seen to lead to much more centrally produced top pairs, and the charge asymmetry has been shown to be a promising observable to distinguish between different new physics models. Our implementation of this new process in POWHEG, called PBZp, is very flexible, as it allows for the simulation of any Z -boson model, and should thus prove to be a useful tool for Z -boson searches in the top-antitop channel at the LHC, in particular for leptophobic models.      where T l denotes the color matrix associated with parton l (T l cb = if clb for gluons, T l ab = t l ab and T l ab = −t l ba for quarks and anti-quarks), s ij = 2p i · p j , and V j (0, 0, ε) = 1 ε 2 , V j (m t , m t , ε) = 1 ε with v ji = 1 − where again s = m 2 t (1 + x) 2 /x and where the double poles are seen to originate only from initial-state massless quarks. The IR poles are given by the Born cross section multiplied by a factor I init + I final .