Associated production of a top pair and a Z boson at the LHC to NNLL accuracy

We study the resummation of soft gluon emission corrections to the production of a top-antitop pair in association with a Z boson at the Large Hadron Collider to next-to-next-to-leading logarithmic accuracy. By means of an in-house parton level Monte Carlo code we evaluate the resummation formula for the total cross section and several differential distributions at a center-of-mass energy of 13 TeV, and we match these calculations to next-to-leading order results.


Introduction
The associated production of a top pair and a Z or W boson are the two processes with the heaviest final states measured to date at the Large Hadron Collider (LHC). The total cross section for these processes was measured during Run I [1,2], and preliminary measurements at a center-of-mass energy of 13 TeV are also available [3,4]. The ttZ production process is particularly interesting because it allows one to study the coupling of the Z boson with the top quark. This measurement further tests the Standard Model (SM) of particle physics and probes several Beyond the SM scenarios that predict changes to this coupling with respect to the SM. In addition, these production processes lead to high multiplicity final states which are background in the search for new heavy states decaying via long chains, such as dark matter candidates.
Given their importance for phenomenological studies, next-to-leading-order (NLO) QCD and electroweak corrections to the associated production of a top pair and a massive vector boson were studied by several groups [5][6][7][8][9][10][11][12]. A full calculation of the QCD corrections to next-to-next-to leading order (NNLO) accuracy would be desirable but it is extremely difficult even with the most up to date techniques for the calculations of higher order corrections. However, the associated production of a top pair and a heavy colorless boson is a multiscale process which is expected to receive potentially large corrections arising from soft gluon emission. The resummation of these effects to next-to-next-to leading logarithmic (NNLL) accuracy can be carried out by exploiting the factorization properties of the partonic cross section in the soft limit (which can be studied with effective field theory methods 1 ) and by subsequently employing renormalization group improved perturbation theory techniques. In the case of the associated production of a top pair and a Higgs boson the resummation formula in the soft emission limit was discussed in [14], and results for the total cross section and several differential distributions at NLO+NNLL accuracy were presented in [15]. Studies of the associated production of a top quark pair and a W boson to NLO+NNLL accuracy can be found in [16], where the resummation was carried out in Mellin moment space as in [15], and in [17], where the resummation was instead carried out in momentum space.
The results of [15] and [16] were obtained by means of an in-house parton level Monte Carlo code for the numerical evaluation of the resummation formula. The output of this code was then matched to complete NLO calculations obtained by employing MadGraph5_aMC@NLO [18] (which we indicate with MG5 aMC in the following). Building on the results of those two papers, in this work we obtain a resummation formula for the associated production of ttZ final state, and we evaluate it to NNLL accuracy by means of dedicated parton level Monte Carlo code. We match our results for the total cross section and differential distributions to NLO calculations in order to obtain predictions at NLO+NNLL accuracy.
The paper is organized as follows: In Section 2 we introduce some basic notation and we briefly summarize the main steps in our calculations. For a more technical discussion of the methods employed in this paper, we refer the reader to the detailed descriptions provided in [14][15][16]. In Section 3 we present predictions at NLO+NNLL accuracy for the total cross section as well as for several differential distributions. Finally, we draw our conclusions in Section 4.

Outline of the Calculation
The associated production of a top quark pair and a Z boson receives contributions from the partonic process where ij ∈ {qq,qq, gg} at lowest order in QCD, and X indicates the unobserved partonic final-state radiation. The two Mandelstam invariants which are relevant for our discussion areŝ The soft or partonic threshold limit is defined as the kinematic region in which z ≡ M 2 /ŝ → 1. In this region, the final state radiation indicated by X in (2.1) can only be soft. The factorization formula for the QCD cross section in the partonic threshold limit is the same as the one derived in [14] for ttH production, up to the straightforward replacement of the Higgs boson with a Z boson: We indicated with s the square of the hadronic center-of-mass energy and we defined τ min = (2m t + m Z ) 2 /s and τ = M 2 /s. The notation adopted for the channel dependent hard functions H, soft functions S, and luminosity functions f f , as well as for the final-state phase-space integration measure, is the same one used in [15,16] and we refer the reader to these papers for more details. Similarly to LO, the only subprocesses to be considered in the soft limit are those labeled by indices ij ∈ {qq,qq, gg}. The hard and soft functions are two-by-two matrices in color space for qq-initiated (quark annihilation) processes, and three-by-three matrices in color space for gg-initiated (gluon fusion) processes. Contributions from other production channels such asqg and qg (collectively referred to as "quark-gluon" or simply "qg" channel in what follows) are subleading in the soft limit. The hard functions satisfy renormalization group equations governed by the channel dependent soft anomalous dimension matrices Γ ij H . These anomalous dimension matrices were derived in [19,20].
In order to carry out the resummation to NNLL accuracy, the hard functions, soft functions, and soft anomalous dimensions must be computed in fixed-order perturbation theory up to NLO in α s . The NLO soft functions and soft anomalous dimensions are the same ones needed in the calculation of ttH and ttW ± to NNLL accuracy and can be found in [14][15][16]. The NLO hard functions are instead process dependent, receive contributions exclusively from the virtual corrections to the tree level amplitudes, and were evaluated by customizing the one-loop provider Openloops [21], which we used in combination with the tensor reduction library Collier [22][23][24][25]. The NLO hard function have been cross-checked numerically by means of a customized version of GoSam [26][27][28][29], used in combination with the reduction provided by Ninja [30][31][32].
In this paper we carry out the resummation in Mellin space, starting from the relation σ(s, m t , m Z ) = 1 2s where f f ij is the Mellin transform of the luminosity functions, and c is the Mellin transform of the product of the hard and soft function (see [15,16] for details). Since the soft limit z → 1 corresponds to the limit N → ∞ in Mellin space, we neglected terms suppressed by powers of 1/N in the integrand of (2.4).
The hard and soft functions included in the hard scattering kernels c in (2.4) can be evaluated in fixed order perturbation theory at scales at which they are free from large logarithms. We indicate these scales with µ h and µ s , respectively. Subsequently, by solving the renormalization group (RG) equations for the hard and soft functions one can evolve the factor c to the factorization scale µ f . Following this procedure one finds Large logarithmic corrections depending on the ratio of the scales µ h and µ s are resummed in the channel-dependent matrix-valued evolution factors U.
The expression for the evolution factors is formally identical to the one found for ttW and ttH production and can be found for example in equation (3.7) of [16]. The l.h.s of (2.5) is formally independent of µ h and µ s . In practice however, one cannot evaluate the hard and soft functions at all orders in perturbation theory; this fact creates a residual dependence on the choice of the scales µ h and µ s in any numerical evaluation of c. The hard and soft functions are free from large logarithms if one chooses µ h ∼ M and µ s ∼ M/N . The choice of a N -dependent value for µ s produces a branch cut for large values of N in the hard scattering kernels c, whose existence is related to the Landau pole in α s . We deal with this issue by choosing the integration path in the complex N plane according to the Minimal Prescription (MP) introduced in [33].
Finally, the parton luminosity functions in Mellin space, which we need in the numerical evaluations, are constructed using techniques described in [34,35].

Numerical Results
The main purpose of this section is to present predictions for the associated production of a top pair and a Z boson to NLO+NNLL accuracy. However, we also analyze sys- ii) Approximate NLO calculations, obtained from the NLO expansion of the NNLL resummation formula. We check that for our choice of scales and input parameters approximate NLO calculations provide a satisfactory approximation to the exact NLO calculation. The approximate NLO formulas obtained by expanding (2.5) account for the single and double powers of ln N as well as N -independent terms but not for terms suppressed by inverse powers of N . Nindependent terms depend on the Mandelstam variables, however we refer to them as "constant" terms in what follows. The approximate NLO formulas are obtained by setting µ h = µ s = µ f in the NNLL version of (2.5). Approximate NLO calculations are carried out with the in-house parton level Monte Carlo code which was developed specifically for this project.
iii) NLO+NLL calculations, which are obtained by matching NLO results with resummed results at NLL accuracy obtained by means of the in-house Monte Carlo code. The results are matched according to the formula The terms in the square brackets, which contribute to NLO and beyond, depend on the scales µ s and µ h in addition to the factorization scale µ f . Of course the dependence on µ s and µ h is formally of NNLL order; by varying these scales and the factorization scale in (3.1) one can estimate the size of the NNLL corrections.
iv) NLO+NNLL predictions are obtained by evaluating the hard scattering kernels in (2.5) to NNLL accuracy with the in-house Monte Carlo code and by matching the results to NLO calculations as follows: The terms in the squared brackets in (3.2) contribute starting from NNLO and represent the NNLL corrections to be added to the NLO result.
v) Approximate NNLO calculations are obtained by the NNLL resummation formula and include all powers of ln N and part of the constant terms from a complete NNLO calculation. The approximate NNLO formulas employed in this paper are constructed as the ones employed in [15,16] for ttW and ttH production. A detailed description of the constant terms which are included in the approximate NNLO formulas can be found in Section 4 of [14]. Approximate NNLO formulas are evaluated with the in-house Monte Carlo code which we developed and they are matched to the NLO calculations as follows where we label the matched result "nNLO" for brevity. By construction nNLO predictions are independent from the hard and soft scales but they do have a residual N 3 LO dependence on µ f . vi) NLO+NNLL expanded to NNLO. Finally we consider a second way of expanding the NNLL resummation formula to NNLO. This approach differs from the approximate NNLO result used above by constant terms, which are formally of N 3 LL accuracy. This approximation is defined by the relation The constant pieces in (3.4) contain explicit dependence on µ h and µ s , in addition to that on µ f . This dependence is formally an effect of N 3 LL order. By comparing the predictions obtained from (3.4) to the corresponding NLO+NNLL calculations one can see the relative weight of terms of N 3 LO and higher in the NLO+NNLL calculations. If in the future a complete NNLO calculation for the ttZ production cross section were to become available, it would be possible to match it to the NNLL resummation formula by using precisely this kind of NNLO expansion of the NNLL resummation, as can be seen by replacing N → NN in all of the superscripts in (3.1).

Scale choices
Since any numerical evaluation of the resummed expression for the hard scattering kernels must be carried out by evaluating the factors in (2.5) up to a certain order in perturbation theory, the resummed kernels c will show a residual dependence on the scales µ s and µ h . In order to follow closely the approach adopted in "direct QCD" calculations [33,37,38], the standard choice which we adopt in this work for the hard and soft scales is µ h,0 = M and µ s,0 = M/N [15,16,39,40].
In addition, both fixed order and resummed calculations depend on the factorization scale µ f , which should be chosen in such a way that the logarithms of the scale ratio µ f /M are not large [41]. It is therefore reasonable to choose a dynamical default value for the scale µ f which is related to M . The dependence of the total ttZ production cross section on the ratio µ f /M at the LHC with √ s = 13 TeV is shown in Figure 1. Each line corresponds to a different perturbative approximation, as indicated in the legend. Figure 1 shows that the NLO, NLO+NLL and NLO+NNLL curves intersect in the vicinity of µ f /M = 0.5 and differ significantly for µ f /M 0.5 and for µ f /M 0.5. Following this observation, the default value that we employ for the factorization scale is µ f,0 = M/2.   The uncertainty related to the choice of the factorization scale in fixed order results is estimated as usual by varying this scale in the range µ f ∈ [µ f,0 /2, 2µ f,0 ]. Resummed results depend also on the hard and soft scales, consequently, the uncertainty of the resummed results is estimated by varying separately all the three scales around their default values in the interval µ i ∈ [µ i,0 /2, 2µ i,0 ] for i ∈ {s, f, h}. The scale uncertainty above (below) the central value of a resummed observable O, which can be the total cross section or the value of the differential cross section in a given bin, is determined as follows. First one evaluates the quantities for i ∈ {s, f, h}. In (3.5) we defined κ i = µ i /µ i,0 , andŌ indicates the observable evaluated at κ i = 1 for all i-s. The scale uncertainty above (below)Ō is then obtained by combining in quadrature ∆O + i (∆O − i ) for i ∈ {s, f, h}.

Total cross section
In this section we analyze the total cross section for the associated production of a top quark pair and a Z boson at the LHC operating at a center-of-mass energy of 13 TeV. The relevant results are collected in Table 2. We first compare the approximate NLO cross section, obtained by expanding the resummation formula to NLO (second row of Table 2) with the complete NLO cross section (fourth row) and the NLO cross section without the contribution of the quark-gluon channel (third row). The difference between the approximate NLO result and the NLO result without the qg channel is due to terms in the quark annihilation and gluon fusion channels which are subleading in the partonic threshold limit. We see that the impact of these terms is around 1%. The difference between these two results is therefore small in spite of the fact that the NLO corrections are large, as can be seen by comparing them with the LO result. However, we see that the approximate NLO result shows a smaller scale uncertainty than the NLO result with the contribution of the qg channel. We conclude that the soft emission corrections provide the bulk of the NLO corrections for this choice of the factorization scale. This motivates us to study the effect of the resummation of these corrections, keeping in mind that by matching the resummed results to NLO calculations we consider both power corrections and the contribution of the qg channel to that order. The NLO+NLL and NLO+NNLL cross sections, shown in the sixth and seventh line of Table 2 are main results of this paper. By looking at the NLO, NLO+NLL, NLO+NNLL results we see that the cross section is progressively increased, but the central value of each prediction falls in the scale uncertainty band of the predictions of lower accuracy.
One can then look at the NNLO expansions of the NNLL resummation formula, which are shown in the last two lines of Table 2. By comparing these results to the NLO+NNLL cross section, one sees that the effect of the resummation corrections beyond NNLO are relatively small. As it was observed in the case of the ttH and ttW processes in [14][15][16], the scale uncertainty affecting the nNLO result is very small compared to the NLO+NNLL scale uncertainty, and most likely underestimates the residual perturbative uncertainty at NNLO.
Experimental collaborations reported measurements of the ttZ total cross section in combination with measurements of the ttW cross section [1][2][3][4]. We conclude this section by comparing our predictions for ttW and ttZ with experimental data. The ttW production cross section was evaluated by running the code developed in [16] with the scale choices and input parameters employed in the present work for ttZ production and described above. The results for ttZ and ttW production cross section at 8 and 13 TeV are summarized in Table 3. In Figure 2 we follow the structure of Figure 12 in [4] in order to compare graphically calculations with the corresponding experimental measurements. The experimental measurements at 8 TeV are taken from [2], while the experimental measurements at 13 TeV are taken from [4]. The green dots and cross-  Table 3. Total cross section for ttZ and ttW production at the LHC with   shaped "error bars" correspond to NLO calculations carried out with µ f,0 = M/2 and their scale uncertainty. The red dots and crosses correspond instead to NLO+NNLL calculations. It is interesting to observe that, while predictions for the ttZ production cross section are in perfect agreement with the measurements at both 8 and 13 TeV, the predictions for the ttW cross section are slightly smaller than measurements for both collider energies. This observation holds for NLO and NLO+NNLL calculations alike. Of course this discrepancy should be taken with a grain of salt, and requires a more detailed discussion with the experimental collaborations. Moreover, we would like to stress that a fully exhaustive comparison between predictions and measurements should also account for the uncertainty associated to the choice of the PDFs and to the value of α s . These two sources of uncertainty are not reflected in the error bars of Figure 2.

Differential distributions
In this section we obtain predictions for four differential distributions which depend on the momenta of the final state massive particles. The distributions are i) the distribution differential with respect to the ttZ invariant mass, M , ii) the distribution differential with respect to the tt invariant mass, M tt , iii) the distribution differential with respect to the transverse momentum of the top quark, p t T , and iv) the distribution differential with respect to the transverse momentum of the Z boson, p Z T . Figure 3 compares the approximate NLO calculations, carried out with our inhouse code, with the complete NLO calculations, carried out with MG5 aMC. We see that the approximate NLO calculations reproduce well the full NLO calculations. The lower part of each panel shows the ratio between the approximate NLO or complete NLO calculations and the central value of the NLO calculation. One can see that the approximate NLO scale uncertainty band is included in the NLO scale uncertainty band. Figure 4 repeats the same analysis but it compares approximate NLO calculations to NLO calculations without the quark-gluon channel contribution. As expected approximate NLO distributions and NLO distributions without the qg channel have the same shape and scale uncertainty bands of similar size. These two figures show that, for this choice of the factorization scale at least, soft emission corrections provide the bulk of the NLO corrections.  and µ h as described above, is smaller than the NLO scale uncertainty band obtained by varying µ f .
Results at NLO+NLL and NLO+NNLL accuracy are compared in Figure 6. The main effect of the NNLL correction with respect to the NLL ones is an increase of the central value of the bins in the tail of the M and M tt distributions. The scale uncertainty bands turn out to be of similar size at NLO+NLL and NLO+NNLL in almost all bins shown.    of the tail of the p t T distribution, where this envelope includes the NLO+NNLL scale uncertainty.

Conclusions
In the present work we carried out the resummation of soft gluon emission corrections to the associated production of a top-antitop quark pair and a Z boson. The resummation was studied in the partonic threshold limit z → 1 and was implemented to NNLL  accuracy. Numerical calculations of the total cross section and differential distributions to NNLL accuracy were carried out by means of an in-house partonic Monte Carlo code which we developed for this work. The output of this code was matched with NLO calculations obtained from MG5 aMC. The final outcome of this work is represented by the NLO+NNLL calculations of the total cross section and differential distributions for the LHC operating at a center-of-mass energy of 13 TeV presented in the previous section. The code can be easily adapted to carry out phenomenological studies which include cuts on the top, antitop and/or Z boson momenta.
With the choice of the factorization scale made in this work, we can conclude that the soft emission corrections to ttZ production evaluated to NNLL accuracy lead to a moderate increase of the total cross section and differential distributions with respect to NLO calculations of the same observables. The residual perturbative uncertainty at NLO+NNLL accuracy, estimated by varying the soft, hard and factorization scales as explained in the text, is smaller than the NLO scale uncertainty, thus making our evaluations of the cross sections and differential distributions in ttZ production the most precise results currently available in the literature.
This work completes a series of papers devoted to the study of the associated production of a top pair and a colorless heavy boson to NLO+NNLL accuracy in the partonic threshold limit. In [16] the associated production of a top pair and a W boson was studied with the methods employed here for ttZ production, while the associated production of a top pair and a Higgs boson at NLO+NNLL accuracy was considered in [15]. In all cases the resummation was carried out in Mellin moment space. The hard and soft scales were chosen in the same way as in the traditional "direct QCD" approach. Codes for the numerical evaluation of the resummation are now available and tested for all of these three processes, and can be further employed in more specific phenomenological studies, according to the interests of the experimental collaborations. Within such interactions with the experimental community, a detailed study of the uncertainty associated with the choice of the PDFs and to the value of α s (M Z ), in the light of a comparison with the new measurements which are expected in the forthcoming months, would be particularly illuminating.
At this stage, it would also be interesting to combine the NLO+NNLL calculations of ttW , ttZ and ttH production with the electroweak corrections for these processes [11,12]. In addition to this, the inclusion of the decays of the heavy particles in the spirit of [42] is also possible. This would allow to put kinematic cuts on the momenta of the detected particles.