Constraining couplings of the top quarks to the Z boson in ttbar+Z production at the LHC

We study top quark pair production in association with a Z boson at the Large Hadron Collider (LHC) and investigate the prospects of measuring the couplings of top quarks to the Z boson. To date these couplings have not been constrained in direct measurements. Such a determination will be possible for the first time at the LHC. Our calculation improves previous coupling studies through the inclusion of next-to-leading order (NLO) QCD corrections in production and decays of all unstable particles. We treat top quarks in the narrow-width approximation and retain all NLO spin correlations. To determine the sensitivity of a coupling measurement we perform a binned log-likelihood ratio test based on normalization and shape information of the angle between the leptons from the Z boson decay. The obtained limits account for statistical uncertainties as well as leading theoretical systematics from residual scale dependence and parton distribution functions. We use current CMS data to place the first direct constraints on the ttbZ couplings. We also consider the upcoming high-energy LHC run and find that with 300 inverse fb of data at an energy of 13 TeV the vector and axial ttbZ couplings can be constrained at the 95% confidence level to C_V=0.24^{+0.39}_{-0.85} and C_A=-0.60^{+0.14}_{-0.18}, where the central values are the Standard Model predictions. This is a reduction of uncertainties by 25% and 42%, respectively, compared to an analysis based on leading-order predictions. We also translate these results into limits on dimension-six operators contributing to the ttbZ interactions beyond the Standard Model.


Introduction
After run I of the Large Hadron Collider (LHC) at √ s = 7 and 8 TeV, we look back on a highly successful research program. Already this first phase of exploring a new energy regime has provided many exciting results: the Higgs boson was discovered [1,2], its quantum numbers and couplings are highly constrained, and many Standard Model (SM) measurements are competitive with previous ones, if not exceeding them. Furthermore, a plethora of searches for signals of new physics have been undertaken, reaching out into the multi-TeV region as well as exploring small deviations of SM parameters. The absence of any spectacular signal of new physics highly constrains many minimal extensions of the SM and, at the same time, opens up new avenues for experimental searches and theoretical model building. These developments represent a remarkably fast progress and demonstrate the potential of the LHC in the years to come. One particularly promising class of SM processes is top quark pair production in association with gauge bosons or a Higgs boson. Due to their relatively high production threshold these processes were not accessible at the Tevatron. In contrast, the high energy and large luminosity of the upcoming LHC runs will produce sufficiently many events to allow detailed studies of these processes. Progress in this direction has already been made with cross section measurements of tt + γ production by ATLAS at 7 TeV [3] and CMS at 8 TeV [4]. First events for the processes tt + Z/W have also been reported in Refs. [5,6]. It is exciting to envision future studies of these processes with direct measurements of the couplings and new sensitivity to physics beyond the Standard Model.
In this paper we focus on the determination of the top quark to Z boson couplings through ttZ production at the LHC. This process is a direct probe of the ttZ interactions, which distinguishes it from other indirect probes such as the LEP measurements of the ρ-parameter [7] and the Z → bb branching ratio [8]. The SM unambiguously predicts the strength of these couplings, and higher order electroweak corrections modify the leading order values only minimally [9]. On the other hand, extensions of the SM which address, e.g. dynamic electroweak symmetry-breaking, typically induce larger deviations. Popular examples are certain variants of Supersymmetry [10,11] or Little Higgs Models [12,13]. More generally, any new fermion which mixes with the third generation quarks might induce deviations to the ttZ SM couplings. Hence also 4th generation quarks [14,15], topcolor models [16] and extra-dimensional extensions of space-time [17] have to be considered. It is therefore important to know to what extent LHC experiments are sensitive to physics beyond the SM in ttZ production. Clearly, this is not only a question of experimental sensitivity but also depends crucially on our theoretical understanding of the production and decay dynamics of the pp → ttZ process.
The ability of LHC experiments to constrain the ttZ couplings was first considered in a series of studies by Baur, Juste, Orr and Rainwater [18,19]. The authors identified suitable observables which are sensitive to vector and axial couplings as well as to the weak electric and magnetic dipole moments. The tri-lepton signature with semi-hadronically decaying top quarks and a leptonically decaying Z boson turns out to provide a good compromise between clean signature and large enough cross section. But even decay modes with a Z boson decaying into neutrinos yield additional sensitivity [19]. These analyses show that sensitivity to the form factor of the vector current is relatively weak and limits can only be placed within a factor of three with respect to the SM value. In contrast, the form factor of the axial current can be pinned down to about 20% accuracy. The authors of Ref. [20] perform a similar analysis using the more modern language of effective operators. This allows them to relate tbW and ttZ couplings in a combined study of single top and ttZ production.
In the context of this work it is important to emphasize that all previous coupling studies were performed at leading-order, and the large residual scale uncertainty was identified [18] as the main obstacle to stronger constraints on the ttZ couplings. It is the aim of this paper to reduce these uncertainties through a NLO QCD calculation for ttZ production and decay into a realistic final state with leptons, jets and missing energy. The hadronic production of ttZ, with stable top quarks and a stable Z boson, was previously calculated at NLO QCD accuracy by Lazopoulos, McElmurry, Melnikov, and Petriello [21], and by Kardos, Papadopoulos, and Trocsanyi [22]. The latter calculation was also interfaced to a parton shower [23], accounting for the decays of the top quarks and Z boson through the spin uncorrelated parton shower approximation. Further hadronization effects were studied in Ref. [24]. Since our coupling analysis relies on studying leptonic opening angles we believe that spin correlations are crucial for a correct interpretation of the results. We therefore account for NLO QCD spin correlations in the decay of top quarks and hadronically decaying W bosons. This includes the full one-loop corrections as well as soft, collinear and wide angle gluon emission off the top quark decay chain. Spin correlations of the leptonically decaying Z boson are included as well. While including all spin correlations, we approximate top quarks and the Z boson as close to on-shell in the narrow-width approximation. This approximation is parametric in Γ/m and its wide range of validity in tt production has been studied in Refs. [25][26][27][28].
It is interesting to note that the ttZ couplings may also be directly probed through single top production in association with a Z boson. Indeed, the inclusive cross section of tZ plus its charge conjugate processtZ is comparable to the inclusive ttZ cross section [29]. It turns out that this process is also the leading background to a ttZ signal, while other backgrounds such as pp → W Zbbjj are almost negligible [18]. However, it is possible to separate ttZ and tZ production by cutting on forward jets and demanding a high jet multiplicity, including two b-tagged jets [29]. We will therefore consider only the process pp → ttZ in this paper, and defer the study of the couplings using tZ (or a combination of both processes) to a later date.
Finally, let us note that a coupling analysis is not the only scenario in which the process pp → ttZ is interesting. The semi-hadronic decay mode of the top quark pairs together with the leptonic Z boson decay is background to several tri-lepton and same-sign lepton searches with additional jets and missing energy. Such signatures can arise from gluino decays in Supersymmetry, in the context of Universal Extra Dimensions, as well as in models with fermionic top quark partners. Furthermore, the invisible decay Z → νν produces a top pair plus a large amount of missing transverse energy, and is therefore an irreducible background to searches for scalar or fermionic top quark partners decaying into top quarks plus dark matter candidates. While we do not address these topics in this paper, it would be interesting to study the effects of NLO corrections when strong selection cuts are applied on this background.

Outline of the calculation
In this section, we briefly discuss the features of our calculation. We consider the tri-lepton signature pp → tt+Z → t(→ νb)t(→ jjb) Z(→ ) which profits from a large cross section due to the hadronic decay of one W boson and the lepton multiplicities from the remaining W and Z bosons. In our results we will sum over all combination of e ± and µ ± in the final state, allowing either t ort to decay leptonically. Application of the narrow-width approximation for top quarks and the Z boson allows us to separate production and decays stage according to where dσ denotes the production cross section and dB X→Y = dΓ X→Y Γ tot X are the partial branching fractions. The use of the narrow width approximation neglects contributions which are parametrically suppressed by O(Γ/m), arising from a largely off-shell top quark or Z boson. Severe selection cuts on final state particles can violate this approximation when distorting the Breit-Wigner line shape of the resonance. In our analysis we aim for a large cross section and only place mild cuts required by experimental detector acceptance.
Hence, we believe the narrow-width is an excellent approximation for our study 1 . We also neglect the contribution from the decay t → W b + Z since the available phase space for on-shell top quarks is tiny and B t→W bZ ≈ 3 × 10 −6 [30][31][32][33].

NLO QCD correction
At leading order, the production of ttZ occurs through the gg and qq partonic channels. At next-to-leading order QCD, these channels receive real and virtual corrections, while real emission corrections open up the partonic channels qg andqg. We also include NLO QCD corrections to the top quark decays and the hadronically decaying W boson; consequently their total widths are included at LO and NLO as well. Eq. (2.1) expanded up to NLO accuracy reads, arises from the α s expansion of the total widths in the denominator. The virtual corrections are evaluated using a numerical OPP realization [34] of D-dimensional generalized unitarity [35][36][37] (for a review, see Ref. [38]). We extended the framework of Ref. [39] to account for color neutral bosons, which requires new tree level recursion relations as well as an extension of the OPP procedure. Soft and collinear singularities in the real emission corrections are regularized using the dipole subtraction scheme of Refs. [40,41], supplemented with a cut-off parameter for the finite dipole phase space [42][43][44][45]. The virtual and real corrections to the top quark decay and hadronic W boson decay are implemented analytically. Soft and collinear singularities in the real emission decay phase space are regularized using subtraction dipoles given in Ref. [46]. We also would like to point out our utilization of parallel computing features. We implemented a version of the Vegas integration algorithm which allows parallelization [47] via the Message-Passing-Interface (MPI) [48]. The observed speed-up in run time scales almost linearly with the number of CPU cores used. This allows us to obtain a full NLO QCD prediction for the total cross section within a few hours on a modern desktop computer with 8 cores.
We perform several checks to ensure the correctness of our calculation. The squared amplitudes for tree level and real emission corrections are checked against MadGraph v.2.49 [49]. The cancellation of poles in D − 4 of dimensional regularization between the virtual corrections and integrated dipoles has been verified for several phase space points. We also checked the finite part of the virtual amplitudes against the automated program GoSam [50] for a few phase space points and find very good agreement. Our framework also allows us to turn the Z boson into an on-shell photon which we used to cross check against the amplitudes of Ref. [46]. At the level of the integrated cross section, we vary the cut-off parameter for the finite dipole phase space by at least one order of magnitude and verify independence from this parameter for the total cross section and kinematic distributions. The interface of production and decay amplitudes is checked by integrating over the full phase space and verifying the factorization into the inclusive cross section for stable top quarks and Z boson times their branching ratios, at NLO QCD. Finally, we compare our full hadronic results with a previous calculation [24] in the literature for stable top quarks and Z boson. For the purposes of this comparison, we take masses of the top quark, m t + m Z /2, we find a leading order cross section of 103.5(1) fb and a next-to-leading order QCD cross section of 137.0(3) fb. This is to be compared with the results of Ref. [24] of 103.5(1) fb and 136.9(1) fb, at leading and next-to-leading order QCD. The cross sections are in excellent agreement within the integration errors. Figure 1 also demonstrates good agreement in shape for the top quark p T distribution between our results and Fig. 1a in Ref. [24].

ttZ couplings
The ttZ interaction Lagrangian in the SM can be written as with the electromagnetic coupling constant e. The vector and axial couplings are where Q t = 2/3 is the top quark electric charge, T 3 t = 1/2, and θ w is the weak mixing angle. The numerical values for the SM couplings are C SM V 0.244 and C SM A −0.601. New physics contributions to the ttZ couplings are most conveniently introduced by higher dimensional operators in the language of effective field theory. A minimal set of dimensionsix operators for top quark production and decay have been categorized in Refs. [53][54][55]. In total there are 91 different operators which can be summarized into 20 different anomalous couplings, if on-shellness and gauge invariance is enforced [54]. For interactions of a Z boson with top quarks only four anomalous couplings, C 1/2,V /A , remain and Eq. (2.3) becomes In this work we will confine ourselves to the study of the above vector and axial couplings C 1,V /A , and neglect the C 2,V /A terms corresponding to the weak magnetic and electric dipole moments of the top quark. Their tree level values vanish in the SM, and C 2,V receives one-loop corrections of O(10 −4 ) [56], while C 2,A receives finite contributions only beyond two-loops [57]. On the more technical side, the tensor structure that multiplies the C 2,V /A couplings introduces the complication of non-renormalizable amplitudes at NLO QCD. While it is straightforward to handle such contributions, our current implementation of the OPP integrand reduction method does not allow tensor ranks larger than N for N -point loop integrals. Such an extension of the OPP reduction algorithm has been outlined in Appendix B of Ref. [58]. We intend to come back to this issue in a separate publication in order to study the phenomenological implication of electroweak dipole moments in ttZ production.
The couplings C 1,V /A can be written in terms of the SM contribution plus deviations due to higher dimensional operators In the above, t R,L are the top quark spinors and φ is the SM Higgs doublet. For further definitions, we refer the reader to Ref. [54].
We now would like to comment on existing constraints on the ttZ couplings. Clearly, these constraints are not obtained directly through on-shell production of a Z boson in association with top quark pairs. Instead, they arise from potential deviations which the higher dimensional operators in Eq. (2.6) introduce to the ρ parameter and the Zbb vertex in the SM. Those parameters are highly constrained through the experimental fits [59] of the ε parameters [60][61][62], The SM predicts their values as ε SM 1 = (5.21 ± 0.08) × 10 −3 and ε SM b = −(6.94 ± 0.15) × 10 −3 [59]. The new physics contributions in Eq. (2.6) introduce the corrections [63] (2.10) The experimentally measured values in Eq. (2.8) can now be used to constrain the operators C and C 33 φu . We will present the numerical results later in Sect. 3.4 together with our results from ttZ production. A further experimental constraint arises from the measurements of the Zb LbL couplings from R b and A b FB at LEP, which are in per-mille level agreement with the SM predictions [8]. This experimental fact together with the SU(2) L symmetry of the SM can be used to relate C . Hence, one of the two operators can be eliminated from Eq. (2.6).

NLO QCD Results
In this section we describe the details of our numerical analysis and the results. We consider the process pp → tt + Z → t(→ νb)t(→ jjb) Z(→ ) and sum over all combinations of leptons e ± , µ ± . We choose the following fixed input parameters Unless otherwise stated, we use MSTW2008 parton distribution functions [64] with α s (M Z ) = 0.13939 and α s (M Z ) = 0.12018 at LO and NLO, which we evolve to the renormalization scale µ ren using 1-loop and 2-loop running, respectively. The LO and NLO scale dependence has been studied in previous works. We do not repeat these studies here, and adopt the central scale µ 0 = m t + M Z /2 for µ = µ ren = µ fact , as suggested in Ref. [21]. Since we include NLO QCD corrections to the top quark decay and the hadronically decaying W boson, we need to include their total widths up to next-to-leading order,  We consider proton-proton collisions at the LHC with a center-of-mass energy of √ s = 13 TeV. To account for detector acceptances and trigger we require Jet are defined by the anti-k T algorithm [65] with R = 0.4. With these input parameters and cuts we find the LO and NLO QCD cross sections, for the central scale µ 0 which is varied by factors of 2 and 1/2. The value in brackets is the integration error on the last digit. The dependence on the unphysical scale is reduced from approximately ±30% at LO to ±13% at NLO QCD. Higher order corrections increase the cross section by 36%, K = σ NLO ttZ σ LO ttZ = 1.36. We also calculate the cross sections without any acceptance cuts and find a significantly lower K = 1.23. This emphasizes the importance of modeling a realistic final state with all unstable particles decayed. The ratio of the cross sections with and without cuts defines the acceptance function A, for which we find The increase of approximately 3% when going from leading to next-to-leading order seems minor. However, the common practice of modeling acceptance effects at LO and multiplying with a K-factor obtained from a NLO calculation with stable particles, underestimates the correct NLO cross section by  a calculation that uses the pdf sets from CTEQ6L1 [51] and CT10 [66] at LO and NLO QCD, respectively. We find These cross sections are about 14% smaller at LO and 7% smaller at NLO QCD compared to the results obtained with MSTW parton distribution functions. At NLO, we find that 4% out of the total difference of 7% is due to the different values of α s (M Z ). The resulting scale uncertainty bands are approximately the same for CTEQ and MSTW pdfs. Hence, the difference due to the two different parton distribution sets is well within the uncertainty estimate from factorization and renormalization scales. Before turning to the ttZ coupling analysis, let us discuss some generic kinematic distributions. Fig. 2 (left) shows the transverse momentum of the two lepton system reconstructing the Z boson. Similar to the total cross sections we observe a strong reduction in unphysical scale dependence over the entire p T spectrum. Scale bands for LO and NLO predictions are comfortably overlapping. From this plot we read off an average transverse momentum of the Z boson of almost 100 GeV with a far-extending kinematic tail, promising approximately 30 events with p Z T ≈ 300 GeV from 300 fb −1 at the 13 TeV LHC. Fig. 2  (right) shows the azimuthal opening angle between the two leptons from the Z boson decay. This observable has been proven to be a good analyzer of the ttZ couplings [18] and we will consider it in the following analysis. Again, we observe strong reduction in scale dependence when going from LO to NLO. The differential K-factor in the lower pane of these plots shows shape changes in the range of 10% due to higher order corrections, which makes the p Z T spectrum harder and decreases the opening angle between the leptons. On the left hand side of Fig. 3 we compare the transverse momentum spectra (at NLO QCD) of the leptons arising from either the top quark decay or the Z boson. The leptons arising from the Z boson are significantly harder than those arising from the top decay. In  z distribution show that the two NLO predictions obtained with MSTW [64] and CTEQ [51,66] pdfs yield consistent results over the entire spectrum. However, as can be seen in the lower pane, the K-factors differ significantly (10% or more) due to very different predictions with LO pdfs (cf. also Eqs. (3.4) and (3.6)).

Coupling extrapolation and statistical analysis
We will now use our calculation to investigate the constraints that can be placed on ttZ couplings, using both existing and anticipated LHC data. To do so, we need to determine how normalization and shapes of differential distributions depend on variations of the couplings. Hence total cross sections and differential distributions need to be calculated for a large grid of C 1,V and C 1,A coupling values. This is simple enough at LO, and while it is still feasible at NLO, it does place a strain on computing resources. As a convenient alternative, we note that ttZ production and decay amplitudes at LO or NLO QCD can be written as with the coefficients M i encoding both the kinematics and all couplings other than the ttZ couplings. The differential cross section is then dependent on six coupling structures, and can be written as Evaluating the cross section for six values of (C 1,V , C 1,A ) allows us to solve for the coefficients s i . These can then be used to extrapolate results for any values of C 1,V and C 1,A . Furthermore, this fitting procedure can not only be done for the total cross section but also bin-by-bin for a given distribution, retaining the effects of spin correlations and selection cuts. As a check of this approach, we have evaluated the cross sections and distributions for a few points in the (C 1,V , C 1,A ) parameter space, both by an explicit calculation and 13 TeV by using the fit for the s i coefficients. Excellent agreement is found in all cases. As an example, we show one comparison in Fig. 4 for the ∆φ + z − z distribution, which we will later use in the coupling analysis. As can be seen in Fig. 4, the overall normalization and the shape are correctly reproduced by the fitting procedure. In this figure and the following, the relative shifts in the couplings are given by In our analysis we focus on the tri-leptonic final state and employ the azimuthal angle between the leptons originating from the Z boson decay to perform our analysis. This angle has been identified as being particularly sensitive to the ttZ couplings in Ref. [18]. We have already discussed the strong reduction in scale uncertainty when going from LO to NLO QCD for this observable. Here, in Fig. 5 (left), we show the effect of NLO QCD corrections on the shape of the normalized ∆φ distribution. Higher order effects tend to shift events from larger to smaller opening angles. In Fig. 5 (right) we show that similar shape changes can arise due to variations of the vector and axial ttZ couplings. This emphasizes the importance of precise predictions, since missing higher order effects might be misinterpreted as deviations from the SM. To illustrate that the ∆φ shape is a useful discriminator for our coupling analysis, we have chosen a value (∆C 1,V , ∆C 1,A ) in Fig. 5 (right) such that the total cross section approximately coincides with the SM ttZ cross section. Hence, a measurement of the rate alone would not reveal the deviations from its Standard Model value.
Let us now outline the basic features of our statistical analysis. We are interested in answering the question: what are the bounds that can be placed on deviations of the ttZ couplings, assuming that the SM is true? Obviously, the answer will depend on the assumed integrated luminosity of the data sample as well as on theoretical and experimental uncertainties. We assume the SM prediction as our null hypothesis H SM (∆C 1,V , ∆C 1,A ) = (0, 0), against which we test alternative hypotheses H alt with (∆C 1,V , ∆C 1,A ) = (0, 0). Alterna-tively, the null hypothesis can be replaced by observed data from LHC experiments to determine the best fit in (∆C 1,V , ∆C 1,A ) parameter space. This also enables the exclusion of parts of the parameter space at a given confidence level. Assuming that the data agree with the SM, the bounds on the ttZ couplings obtained in this work should approximate those obtained from real data. We construct two likelihood functions L SM and L alt which allow us to define a test statistic Λ = log (L SM /L alt ). We then generate two event samples for a fixed integrated luminosity assuming that either H SM or H alt is true. The test statistic Λ can be evaluated for these two event samples, and repeating this evaluation in a large number of pseudo experiments provides the probability distributions P (Λ|H SM ) and P (Λ|H alt ). The overlap of these two probability distributions can be used to define the type-I error for rejecting H SM in favor of H alt , even though H SM is true. This error can finally be translated into the more familiar confidence level in terms of standard deviations.
Let us now describe the procedure outlined above more precisely and illustrate how differential distributions at NLO QCD can be used. We closely follow typical likelihoodbased analyses as described for example in Ref. [67], based on the original procedure by Feldman and Cousins [68]. The starting point is the binned likelihood function The log-likelihood ratio is guaranteed to be the optimal test statistic through the Neyman-Pearson lemma [69]. It can now be evaluated with ν i being the events from the ∆φ ll histogram, and n obs being the Poisson distributed events around ν i , for either the SM or the alternative hypothesis. Repeating this procedure for a large number of pseudoexperiments yields the two probability distributions of Λ( n SM ) and Λ( n alt ). An example of two such probability distributions, P (Λ|H SM/alt ), is shown in Fig. 6 for LO (left) and NLO QCD (right). These two distributions can be used to define a confidence level for excluding the alternative hypothesis. For a given valueΛ, the probability of accepting H alt even though H SM is correct (type-I error) is Similarly, the probability of accepting H SM even though H alt is correct (type-II error) is given as (3.14) We defineΛ such that α = β, i.e. there is equal chance of incorrectly rejecting one hypothesis in favor of the other. The value α(Λ) is then a measure of statistical discrimination between the two hypotheses. It can be translated into the more familiar number of standard deviations by σ = √ 2 erf −1 (1 − α), (3.15) where erf −1 is the inverse error function. The above discussion made no mention of systematic uncertainties. In this work we would like to include the leading theoretical uncertainties from unphysical scale dependence and errors associated with parton distribution functions. For simplicity we neglect experimental systematics such as efficiencies or momentum smearing effects. Note however that we include realistic detector acceptances through the cuts in Eq. (3.3). Statistical fluctuations are obviously included in our analysis through the Poisson distribution in Eq. (3.10). Following Ref. [70], we include the theoretical uncertainties through nuisance parameters by including multiplicative factors in the log-likelihood function. The inclusion of nuisance parameters removes the Neyman-Pearson guarantee that the log-likelihood ratio is the optimal test statistic. Nevertheless, one still expects the test to be approximately optimal as long as the nuisance parameters are reasonably constrained. We include the theoretical uncertainties by modifying the likelihood function in Eq. (3.10) according to We choose G to be a constant normalized function with supportν H min/max,i =ν H i (1 ± ∆ theor. unc. ), where θ(x) is the unit step function. The value ofν H i is determined by the assumed luminosity times the cross section in the ith bin for the central scale choice µ 0 . To be most conservative, we choose the cross section within the uncertainty band for each hypothesis such that their difference is minimized. The value in each bin is then rescaled accordingly. This treatment results in a larger overlap between the two likelihood distributions, and consequently a larger α value and less discriminatory power between the two hypotheses. This feature is clearly visible in Fig. 6, when comparing the solid with the dotted curves. Contrasting the LO results in Fig. 6 (left) with the NLO result (right) shows that the lower uncertainty associated with the NLO prediction allows for significantly better statistical discrimination between the hypotheses. Additionally, the increase of the NLO cross section due to the K-factor of approximately 1.4 leads to a larger number of expected events and therefore to smaller statistical uncertainties.

ttZ coupling constraints from current and future LHC data
We now apply the analysis outlined in the previous section to study coupling constraints from current and future LHC data. Figure 7 shows the relative shift of the ttZ cross section for non-SM couplings with respect to the SM cross section for a wide range of vector and axial couplings. The grid of 3200 NLO QCD cross sections is generated with the fit described in Eq. (3.8), at low computational cost, and accounts for selection cuts of Eq. (3.3). We find that within the given range the cross sections vary by about ±50% away from the SM value due to shifts of vector and axial couplings. The remaining scale uncertainty at NLO QCD was found to be ±13%, which roughly corresponds to the area enclosed by the dotted line in Fig. 7. Hence, for all coupling values within this band, a rate measurement alone is not sensitive to any deviation. This is true for a large range of couplings far off the SM value, e.g. (∆C 1,V , ∆C 1,A ) = (1.7, −0.3). We will later see that adding shape information from kinematic distributions will improve this situation and lead to a more powerful discrimination. It is clearly noticeable in Fig. 7 that cross sections are symmetric around the axis ∆C 1,V = −1. This feature can be easily understood from the fact that the LO cross section is dominantly proportional to C 2 1,V + C 2 1,A and ∆C 1,V = −1 corresponds to the point C 1,V = 0. We expect to see a similar symmetry around ∆C 1,A = −1, however the sign of the axial coupling is already constrained from LEP measurements of the Zb LbL interaction when SU(2) L symmetry is invoked, and consequently we do not show the results for a negative value of C 1,A here.
As a side remark we would like to briefly note that we also studied the effects of ttZ coupling shifts on the top quark forward-backward asymmetry (A tt FB ) at the Tevatron. At leading order, we considered the parity-violating process qq → Z/γ * → tt with coupling variations as shown in Fig. 7. We find that within this range the forward-backward asymmetry is not significantly enhanced. Hence, any discrepancy between theory and experiment for A tt FB cannot be explained by deviations of the ttZ couplings as assumed in this paper.
We now use current LHC data to obtain direct constraints on vector and axial couplings. The production of ttZ has been observed at the √ s = 7 TeV run at the LHC, with CMS observing nine events [6], and ATLAS observing one event with more stringent selection criteria [5]. This enables ATLAS to place an upper bound on the ttZ cross section, while CMS is able to determine σ CMS ttZ = 0.28 +0.14 −0.11 (stat.) +0.06 −0.03 (sys.) pb. Clearly, error bars from this very first measurement are large, nevertheless it constitutes a 3.3 standard deviation from the background-only hypothesis. The measured total cross section is also consistent with the NLO QCD predictions of 0.137 pb ±11% by Ref. [23] (CTEQ pdfs) or with 0.148 pb ±11% from our calculation (MSTW pdfs). In spite of the low number of events and correspondingly high statistical error, it is instructive to use this measured cross section to place bounds on the ttZ couplings. This constitutes the first direct constraints on these couplings. We perform a log-likelihood ratio analysis, as described in Section 3.2. Our null hypothesis is derived from the experimental cross section σ CMS ttZ = 0.28 pb, from which we predict 24 events in the tri-lepton channel from 5 fb −1 of data if no acceptance cuts are applied. We are adopting this number as our reference point against which we compare the alternative hypothesis, which is a given point in (C 1,V , C 1,A ) parameter space, and for which we calculate a predicted number of events in a 5 fb −1 dataset. To account for uncertainties in the theoretical prediction we include a uniform distribution spanned by the theoretical error band, for which we choose ∆ theor. unc. = 40% at LO and ∆ theor. unc. = 15% at NLO QCD. For experimental systematics we need to introduce a Gaussian distributed probability. Hence, the function G(..|..) in Eq. (3.16) has to be modified and becomes whereν is the mean value of the experimental measurement andσ = σ exp. syst. is the systematic experimental uncertainty. In this way the mean value ν in the likelihood function is normal distributed during the generation of pseudo experiments with a standard deviation σ exp. syst. . We adopt the experimental systematic error of σ exp. syst. = ±20% quoted in Ref. [6]. Again, we use Eq.  Figure 8, for our LO (left) and NLO (right) calculations. The color code shows the significance with which an alternative coupling hypothesis can be excluded with respect to to the experimental data. In the plane of relative deviations of vector and axial couplings, the point (∆C 1,V , ∆C 1,A ) = (0.0, 0.0) corresponds to the SM value. Clearly, this point is fully consistent with the experimental measurement. By comparing left and right plots we notice the stronger constraints when NLO input is used. The constraints from the data using a LO calculation are −11 ∆C 1,V 10 and −4 ∆C 1,A 2 at the 95% C.L., while they improve at NLO to −8 ∆C 1,V 7 and −3 ∆C 1,A 1. Of course, these limits are extremely loose, and furthermore should be interpreted with care since very few events have been observed by the experiments so far. Only a larger data set and detailed analysis of backgrounds and detector effects will provide more reliable constraints on the ttZ couplings. We nevertheless believe that these results are interesting to consider, especially when put into context with limits obtained from the future high-energy LHC.
We turn now to this anticipated run of the LHC at √ s = 13 TeV. Using the interpola- in Fig. 9 show the significance with which non-SM ttZ couplings can be separated from the SM hypothesis, assuming that the SM hypothesis is true. Clearly, this significance is a function of the accumulated luminosity and the associated uncertainties at the given order in perturbation theory. We therefore present six scenarios for luminosities of 30 fb −1 , 300 fb −1 , and 3000 fb −1 at the 13 TeV LHC with theory input at leading and next-to-leading order in QCD. The LO pdf uncertainties are slightly smaller than at √ s = 7 TeV, allowing us to use an overall scale uncertainty of 30% at LO and 15% at NLO. The couplings outside the light-blue area in Fig. 9 roughly correspond to the ones that can be excluded at 68% confidence level (C.L.), whereas couplings outside the orange colored boundary can be excluded at 95% C.L. From comparing the first, second and third row of plots in Fig. 9 it is immediately apparent that increasing the luminosity drastically improves the limits. By comparing plots in the left versus the right column we also see that the bounds at NLO QCD are far stronger. This is a result of the reduced scale uncertainty and the larger cross section due to a positive perturbative correction. Numerically, one finds that with 300 fb −1 and LO input, ∆C 1,V is constrained between −4.0 < ∆C 1,V < 2.8 and ∆C 1,A is constrained between −0.36 < ∆C 1,A < 0.54, at the 95% C.L. 2 The limits improve with NLO QCD predictions to −3.6 < ∆C 1,V < 1.6 and −0.24 < ∆C 1,A < 0.30. In terms of absolute values, these intervals correspond to C V = 0.24 +0.39 −0.85 and C A = −0.60 +0.14 −0.18 at NLO QCD. This is a reduction by 25% and 42%, respectively, compared to results obtained at leading order. A noticeable feature in the exclusion limits at NLO, Fig. 9 (right), is the turnover from a downwards bend shape to an upwards bend shape when going from 30 fb −1 to 3000 fb −1 . This feature is a complicated effect of our uncertainty treatment and a transition from normalization to shape sensitivity in the exclusion. In the upper right plot the exclusion region roughly follows the shape already observed in Fig. 7. This can be understood from the fact that with a small event sample the exclusion limit is dominated by normalization differences, whereas different shapes have vanishing exclusion power. Using a larger event sample, shape sensitivity increases and allows us to exclude regions where normalization differences are small but shapes differ significantly. This is the case in the lower right plot of Fig. 9, where the previous downwards bend is safely excluded.

Limits on dimension-six operators
Having presented our main results in Fig. 9, we can use the obtained limits to put constraints on possible effects of physics beyond the SM. The relevant dimension-six operators have been presented in Sect. 2.2. This is also where we pointed out that the excellent agreement between experiment and prediction for the Zb LbL couplings can be used (together with SU(2) L symmetry of the SM) to relate C . In the following we will make use of this fact and eliminate C (1,33) φq from our analysis 3 . Hence we are left with only two dimension-six operators, C (3,33) φq and C 33 φu . We begin by using the total cross section measurement of CMS at 7 TeV (see Sect. 3.3). The limits on the total ttZ cross section as a function of ∆C 1,V and ∆C 1,A directly translate into limits on the two operators. Diagonalizing the dependence in Eq. (2.6), we find at next-to-leading order This result should be considered with care given the low number of events observed in current experiments. More reliable and stringent limits can only be obtained once more 2 We checked that these LO limits roughly agree with the ones quoted in Ref. [18] for the 14 TeV LHC. 3 It should be noted however that models exist which give vanishing corrections to Zbb for finite C . One example is given in Ref. [71] with vector-like quarks. In such case, our limits remain valid upon the replacement C data is accumulated. To estimate how limits will improve in such a case, we use the results presented in Fig. 9 for the luminosities 30, 300, and 3000 fb −1 . Recall that these results are not only based on the total cross section but also on the shapes of the ∆φ + Due to the weak correlation between vector and axial coupling limits observed in Fig. 9 and because of Eq. (2.6), the limits on C (3,33) φq and C 33 φu are strongly correlated. Hence, the results in Eqs. (3.20)-(3.21) are very conservative. A more appropriate graphical representation of these limits is given in Fig. 10 for our NLO results. We also include the indirect constraints from electroweak precision observables ε 1 and ε b [60][61][62], updated to account for M H = 125.6 GeV in Ref. [59]. All effective operators outside the colored ellipse in Fig. 10 can be excluded at the 95% confidence level. We observe that the limits from ttZ production at the LHC are well-aligned with the precision limit from ε 1 . This can be understood from the fact that ε 1 is directly proportional to the SM ρ-parameter which receives sensitivity from the Z boson self energy with a top quark loop insertion. The constraint from ε b on the other hand arises from the measurement of Z → bb and SU(2) L symmetry of the SM. Hence it leaves C 33 φu mostly unconstrained since this operator contributes to the right handed current only. Altogether, electroweak precision observables put very strong constraints on the ttZ coupling; however, these only arise through indirect sensitivity. Only the analysis of pp → ttZ at the LHC will allow the placing of direct limits for the first time. It would be an interesting continuation of this work to superimpose Fig. 10 with constraints obtained from single tZ +tZ production at the LHC.

Conclusion
In this article we studied top quark pair production in association with a Z boson. Due to its relatively high production threshold and penalties from small branching fractions, this process was never observed at the Tevatron. Even at the 7 and 8 TeV run of the LHC only a few candidate events were collected. As a consequence there is no direct measurement of the top quark to Z boson couplings to this date. This situation will change once the high energy LHC delivers its first tens fb −1 of data. We therefore study the process pp → ttZ in the tri-lepton final state, which provides the best compromise between clean signature and large enough cross section. The central question that we try to answer is by how much limits on ttZ couplings improve once NLO QCD predictions are used. A particularly sensitive observable for such a study is the opening angle between the two leptons from the Z boson decay. We perform the analysis with a binned log-likelihood ratio test which proves advantageous for several reasons. Firstly, the use of likelihood functions guarantees reliable results even for low number of events when, for example, a simple χ 2 analysis would fail. Secondly, non-Gaussian systematic errors such as theoretical scale uncertainties can be straightforwardly implemented in the likelihood ratio test. In addition it turns out to be relatively easy to implement this procedure at NLO QCD. We begin by performing this analysis on the inclusive cross section reported by CMS, which allows us to place the first direct constraints on the ttZ couplings at the LHC. We proceed to an analysis at the √ s = 13 TeV LHC run. Assuming a residual theoretical uncertainty of 15% at NLO we find that with 300 fb −1 of data the vector and axial couplings can be constrained to C V = 0.24 +0.39 −0.85 and C A = −0.60 +0.14 −0.18 at the 95% C.L. This is a significant improvement compared to an analysis at leading order. Even a first determination with only 30 fb −1 of data might be possible if NLO input is used, yielding limits which are about two times weaker. We also translate our constraints on vector and axial couplings into limits on dimension-six operators contributing to the ttZ couplings beyond the SM. The viable region for these operators can be significantly reduced with measurements of pp → ttZ and O(100)fb −1 of data. This allows us to contrast high precision indirect limits from electroweak observables with a direct determination from the LHC.
Finally, we note that effects of New Physics can modify the ttZ coupling beyond vector and axial currents through q 2 -dependent higher dimensional operators. Those couplings typically introduce non-renormalizable interactions and require an extension of our oneloop integrand reduction method. This is an interesting subject for a continuation of this work. Another interesting future topic is the study of sensitivity at a 100 TeV pp collider or at an e + e − machine. At any rate, we look forward to the first analysis of the tt+Z, W, γ, H processes and the subsequent precision studies of top quark phenomenology.