Zγ production at NNLO including anomalous couplings

In this paper we present a next-to-next-to-leading order (NNLO) QCD calculation of the processes pp → l+l−γ and pp→νν¯γ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ pp\to \nu \overline{\nu}\gamma $$\end{document} that we have implemented in MCFM. Our calculation includes QCD corrections at NNLO both for the Standard Model (SM) and additionally in the presence of Zγγ and ZZγ anomalous couplings. We compare our implementation, obtained using the jettiness slicing approach, with a previous SM calculation and find broad agreement. Focusing on the sensitivity of our results to the slicing parameter, we show that using our setup we are able to compute NNLO cross sections with numerical uncertainties of about 0.1%, which is small compared to residual scale uncertainties of a few percent. We study potential improvements using two different jettiness definitions and the inclusion of power corrections. At s=13\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \sqrt{s}=13 $$\end{document} TeV we present phenomenological results and consider Zγ as a background to H → Zγ production. We find that, with typical cuts, the inclusion of NNLO corrections represents a small effect and loosens the extraction of limits on anomalous couplings by about 10%.


Introduction
The power of the Large Hadron Collider (LHC) is its ability both to search for new phenomena at the highest energies and to accumulate a wealth of precise data on a broad range of final states at lower energies. The former is enabled by the center of mass energy of 13 TeV while the latter relies on both the immense amounts of data collected and the technical capabilities of the detector, as well as the ingenuity of the experimental analyses. One of the foremost targets for high-precision studies is the electroweak sector that can be explored, for instance, through the multiple vector boson production channels. Among these, the production rate for a Z-boson and a photon is one of the highest and, at least for leptonic decays of the Z-boson, provides a final state whose particles can all be measured with excellent precision [1][2][3][4][5][6][7][8][9][10][11][12][13]. It is important to perform a careful comparison with the theoretical prediction of the SM since any deviation from it could indicate the presence of anomalous Zγγ or ZZγ couplings [14]. Beyond such tests, Zγ production provides important backgrounds to many other searches. For instance, it is a crucial ingredient in background estimations for the rare Z(→ ll)γ Higgs decay [12]. The kinematics of signal and background are very challenging in this case and require precise predictions to ensure an adequate understanding of both processes [15,16]. The Zγ process also represents a leading background in searches for heavy resonances that decay into a Z-boson and a photon [7,12,13], for example new singlet scalars with loop-induced decays [17]. Beyond that, the channel in which the Z boson decays to neutrinos gives rise to events containing a photon and missing energy. This is therefore of great interest to searches for dark matter and more general BSM studies with similar signatures [18][19][20][21].

JHEP11(2017)150
In order to provide the strongest point of comparison with the SM it is therefore essential that the theoretical prediction for Zγ production be available at the highest possible accuracy. The first steps towards providing a reliable calculation of this cross section, and related observables, were made by extending the initial LO calculations [22,23] to NLO QCD [24]. Further refinements were provided later on, to account for lepton decays in the narrow width approximation and to account for the effects of anomalous couplings [25,26]. The calculation of one-loop helicity amplitudes for this process [27] allowed the inclusion of such effects in a flexible Monte Carlo code [28] and in the generalpurpose code MCFM [29]. The one-loop gluon fusion contribution, gg → Zγ formally enters at NNLO QCD but can also be computed separately [30][31][32] and is important for LHC phenomenology [33]. The effects of NLO electroweak (EW) corrections have also been computed for an on-shell Z-boson [34][35][36] and have been recently included in a full NLO QCD+EW prediction for this process [37]. Finally, the first predictions for Zγ production through NNLO QCD -providing robust theoretical precision for this cross section and related quantities -have been computed recently [38,39].
In this publication we provide an independent calculation of Zγ production at NNLO in QCD. Given the importance of this process, an independent verification such as the one we provide here is invaluable. 1 We extend the previously available results to also include the case of anomalous Zγγ and ZZγ couplings at NNLO. This allows them to be probed with greater confidence and enables more accurate limits to be placed on their possible values. Note that we use the label Zγ here and in the following as a shorthand for calculations that represent the final states with charged leptons l + l − γ and neutrinos ννγ. In particular, in the former case we include the effects of photon radiation in the Z boson decay and also virtual photon contributions. To avoid infrared singularities arising from the emission of photons from partons we use the Frixione smooth cone isolation criterion [43], which we define and discuss in the following section.
Our calculation is based on matrix elements that, for numerical efficiency, have been computed analytically and that are combined to provide a full NNLO calculation using the jettiness phase space slicing method [44,45]. This type of approach to the computation of NNLO corrections is in common with the previous calculation [38,39], that employed a slicing method known as q T -subtraction [46]. These techniques rely on the choice of a slicing parameter that is sufficiently small that either the associated power corrections are negligible or they can be extrapolated away. We perform a detailed check of the jettiness slicing procedure, which is especially important since the numerical uncertainties achieved with slicing methods can be of a similar size as the scale uncertainty observed at NNLO. Focusing on the technical aspects of our calculation is important to validate our results so that they can be made publicly available in the next release of MCFM.
In section 2 we describe the ingredients of our calculation and how we obtain the results presented in the rest of the paper. Section 3 provides a detailed comparison with the results previously obtained in ref. [39]. With our method validated, we proceed to a discussion of JHEP11(2017)150 SM phenomenology in section 4, discussing both Z-boson decays into neutrinos and into charged leptons as an application for the Higgs production background. In section 5 we assess the impact of NNLO corrections on the extraction of anomalous coupling limits. Finally, our conclusions and outlook are given in section 6.

Calculation and setup
Computing Zγ to NNLO accuracy requires the calculation and assembly of numerous contributions that we list here. The existing 0-jettiness slicing implementation in MCFM [47], reviewed briefly below, is used to easily combine all the necessary amplitudes to produce a full NNLO code.
We have taken the Zγ 0/1/2-loop amplitudes from ref. [48], using the function tdhpl [49] to evaluate the two dimensional harmonic polylogarithms through which they are expressed. The double radiation tree amplitudes and the one-loop single radiation amplitudes are based on the previous implementations in MCFM [50] (see also [51]), but are completely rewritten for readability and run-time performance. 2 For the charged lepton decay channel the radiation of a photon from the leptons is possible, which requires the evaluation of the two-loop quark formfactor. We have taken it from ref. [52], where it is given up to three loops, see also [53,54]. Note that the two-loop result has previously already been published in ref. [55].
We have checked that our tree level and one-loop matrix elements agree with results obtained from OpenLoops [56]. Additionally, we found agreement with MadGraph 4.2.7 [57] for our one-loop real emission matrix elements. We work consistently in the five flavor scheme with a zero bottom quark mass.
The current release version 8.0 of MCFM includes Zγ amplitudes with anomalous couplings at NLO [28]. At NNLO one requires the double virtual corrections, easily obtained from the two-loop quark form factor. We additionally implemented the double real radiation anomalous coupling tree amplitudes using the vertices from ref. [58] (following the convention with an additional factor of i as in ref. [59]). Specifically we contracted the vertices with the four parton plus Z boson off-shell current [60,61] 3 to obtain the complete amplitudes. To obtain the real emission one-loop amplitudes we contracted the anomalous coupling vertices with the V → 3j loop amplitude current from ref. [63] and performed the analytic continuations to 2 → 2 kinematics as outlined in ref. [64].
Jettiness subtractions. Since the assembly of our calculation relies on the method of jettiness subtractions, which has been implemented in MCFM-8.0 for color singlet final states [47], we give a brief overview of this method here. For a more complete review we refer the reader to the referenced literature. The N -jettiness subtraction formalism is a phase space slicing scheme for combining infrared singular matrix elements entering higher order perturbative QCD calculations that contain N colored final state partons at Born 2 We have fixed a bug affecting the radiation in decay amplitudes in the existing Zγj implementation in MCFM [50].
3 The hadronic Z current S(+; − : +; −)Ȧ B given in the appendix of ref. [60] for the qqgg case has a typo. The current given in ref. [62] agrees with our own calculation.

JHEP11(2017)150
level [44,45]. Using the event shape variable N -jettiness [65] , τ N , the formalism allows for handling processes with initial-and final-state infrared singularities and any number of colored final state particles. The calculation of the cross section proceeds, not fundamentally different from classic NLO slicing schemes, through a separation of the infrared and collinear singular phase space using a jettiness cutoff parameter τ N cut . The below-cut part is predicted by soft collinear effective theory (SCET) in terms of a factorization theorem using the hard scattering matrix element and soft and beam functions. 4 Here, for Zγ production, we can restrict ourselves to the case of N = 0 colored final state partons at Born level and refer to τ ≡ τ N =0 in the following. For a parton scattering event q a + q b → X + p 1 + . . . + p M with parton momenta q a , q b , p 1 , . . . , p M and color singlets X, the event shape variable 0-jettiness is defined by where q i define the beam and 0-jettiness axes and Q i are in principle arbitrary normalization factors. For example in MCFM we set Q i = 2E i where E i is the energy of an incident parton, such that τ itself is independent of E i and simply depends on the normalized beam axes and p k . 5 For a NNLO cross section, by demanding τ < τ cut the doubly unresolved region can be isolated and the real emission matrix elements in that region can be integrated over analytically in the soft/collinear approximation and finally added to the virtual contributions. The part τ > τ cut corresponds to the singly unresolved limit of the NLO process with an additional jet. It can be handled by any NLO subtraction scheme, preferably using a fully local NLO subtraction scheme like Catani-Seymour subtractions [67,68]. The sum of the two contributions yields a result that, in the limit that τ cut → 0, returns the complete NNLO cross section.
Unless specified otherwise for comparison reasons, we use the boosted jettiness definition [71] which is defined in the color singlet center of mass frame, rather than the hadronic 4 We note that factorization theorems in SCET with respect to other variables also constitute a powerful method for performing higher order calculations; for example top-quark decay has been calculated in such a way at NNLO [66]. 5 In ref. [45] this choice of Qi is referred to as a "geometric measure".

JHEP11(2017)150
definition that is used in MCFM-8.0 and defined in the hadronic center of mass frame. See refs. [45,65,[72][73][74] for discussion. To avoid infrared singularities arising from the emission of photons from partons we use the Frixione smooth cone isolation criterion [43]. In this method one defines a cone of radius R = (∆η) 2 + (∆φ) 2 around the photon where ∆η and ∆φ are the pseudorapidity and azimuthal angle difference between the photon and any parton. The total partonic transverse energy inside a cone with radius R is then required to be smaller than for all cones R < R 0 , where E γ is the transverse photon energy, and , R 0 and n are parameters. When comparing with the results of ref. [39] we adopt their choice of the photon isolation parameters, = 0.5, n = 1 and R 0 = 0.4. For the other results in our study we use the isolation parameters = 0.1, n = 2 and R 0 = 0.4 that for isolated prompt photon production at NLO provide results similar to those obtained using fragmentation functions and experimental cuts [41,75].
Slicing cut-off parameter extrapolation. The leading behavior of the NLO and NNLO cross section coefficients ∆σ NLO , ∆σ NNLO for τ cut → 0 is given by [45], It is these forms that we will use for fitting. We will treat the scale Q as well as the coefficients c i as nuisance parameters that are fitted with the data but unused. See for example ref. [71] for a more detailed discussion of the higher power corrections and their fitting. In all cases we weight the data points in the non-linear fitting procedure with the inverse square of their Monte Carlo integration numerical uncertainty. Our quoted uncertainties for results that have been extrapolated in the limit τ cut → 0, ∆σ NLO and ∆σ NNLO , correspond to a 95% confidence level. Including higher power corrections is equivalent to a theoretical uncertainty in the fit model. To be of use for the prediction, these additional nuisance parameters would require much smaller statistical uncertainties than we demand for the phenomenology in this study. They would become important when their nominal contribution becomes comparable with the statistical uncertainties in our predictions [71]. In any case we make sure that within our uncertainties the results are unchanged by excluding data points with our largest and smallest considered τ cut values. We discuss the extrapolation in more detail in section 3 using actual examples.

Improvements to jettiness slicing
Before we present our main results in the following sections, we first investigate potential improvements to our implementation of jettiness slicing from two sources. Since these improvements are best studied in light of a cutoff independent calculation, we compare NLO
results to those obtained by Catani-Seymour subtractions. The first improvement may be obtained by including the effects of subleading power corrections [71,76] for Drell-Yan type color-singlet production. The second improvement results from a boosted definition of the jettiness slicing variable. We consider the case of a Z-boson decaying into charged leptons at √ s = 8 TeV using the cuts shown in table 1. Figure 1 (left) shows the approach of the jettiness calculation of the NLO coefficient to the result obtained using Catani-Seymour (C.-S.) subtraction [67,68], both with and without the inclusion of the dominant subleading power corrections for the hadronic jettiness definition used in MCFM-8.0. We see that the inclusion of the subleading power corrections has almost no effect. We surmise that this is due to the fact that, in contrast to the processes studied in refs. [71,76], Zγ production contains t-channel diagrams at Born level. To investigate this issue further, we repeat this analysis under a set of cuts that decreases the importance of such diagrams relative to Drell-Yan type s-channel contributions in which the photon is radiated from a charged lepton. This set of cuts corresponds to those specified in table 1, with the replacements: The results of this study are shown in figure 1 (right). Without power corrections the approach to the correct result is similar to the approach when using our standard cuts that is shown in the left panel. However, the inclusion of power corrections now makes a much more marked effect and significantly increases the value of τ cut that might be considered asymptotic. We therefore conclude that the effect of subleading power corrections is process-specific and that the inclusion of the results of refs. [71,76] provides no significant benefit under the cuts considered here. Therefore, and for consistency, we do not include them in any of our subsequent results.
We now compare the performance of the jettiness slicing procedure with two different definitions of τ cut . These correspond to the version in the hadronic center of mass frame employed for other color-singlet processes in MCFM-8.0 and the boosted one in the color singlet center of mass frame. The results of this comparison are also shown in figure 1. We see that the improvement when using the boosted definition is relatively small for the standard cuts, due to the fact that our set of cuts demand the production of a Zγ system that has a rather high virtuality and is quite central. Despite this, the small improvement in the approach to the asymptotic value is the result of a trivial change in the code. We therefore adopt the boosted definition for all the studies presented in this paper.  For a more refined comparison it is instructive to consider not just a total cross section but a differential prediction. Figure 2 shows the pseudorapidity distributions of the NLO coefficients computed using our loose cuts with Catani-Seymour subtraction and using jettiness slicing with τ cut = 1.6 GeV for both boosted and hadronic definitions. For the integrated cross section, as can be seen in figure 1, the difference between the boosted and hadronic definition is about 2.5%. Nevertheless, the difference for rapidities greater than about one already exceeds 2.5% and grows to about 10% for rapidities larger than two. The difference of the integrated cross section between using C.-S. and boosted jettiness is about two percent that manifests as a relatively constant difference throughout the whole rapidity range. This clearly shows that when computing arbitrary differential distributions one cannot generally rely on the smallness of the residual τ cut dependence of the integrated cross section, but one has to consider it differentially as well. For the differential distributions considered in our study we are not affected by this problem, as we show in section 4. This is partly due to directly using the boosted definition, and additionally to just using asymptotically small enough τ cut values where the numerical integration uncertainty dominates over the residual τ cut error.

Validation
The process of Zγ production has been considered previously at NNLO in refs. [38,39] in the q T subtraction formalism [46]. For comparison and mutual validation we check some results of ref. [39] in this section. The 8 TeV cuts used in this comparison have already been summarized in table 1 for the charged lepton decay channel; the counterparts for the neutrino decay channel are shown in table 2.
Specifically we will compare our results with those provided in tables 3 and 5 of ref. [39], for the 8 TeV jet inclusive case. Their stated conservatively-estimated NNLO numerical uncertainty is between 0.5% and 0.6%. Note that the authors do not give a technical explanation of how their extrapolation is performed and the uncertainty is estimated, but only state that it is a combination of statistical Monte Carlo integration error and an extrapolation procedure. As such, our results can evaluate whether their uncertainty estimation is robust.
In order to minimize the uncertainty resulting from the jettiness slicing method we compute the NNLO cross section coefficient separately from the NLO cross section, which is obtained by Catani-Seymour dipole subtractions. From the quoted NNLO results in ref. [39] we subtract our NLO results obtained with NNLO PDFs in order to translate them into an expectation for the NNLO coefficient. Since our NLO results mutually agree within 0.2%, we consider such a procedure for extracting the NNLO coefficients to be reasonable.
We begin with a comparison in the neutrino decay channel, where the authors of ref. [39] report cross sections of 42.33 fb, 70.98 fb and 80.82 fb at LO, NLO and NNLO, respectively. Their stated uncertainty estimate is 0.5% for the NNLO result from inte-

JHEP11(2017)150
Neutrinos  gration error and finite q T cut, which translates to 0.4 fb in absolute terms. We assume that their uncertainty for the LO and NLO results corresponds to the stated number of significant figures, translating to an uncertainty of 0.01 fb. Our LO cross section of 42.35 fb is in good agreement to better than half a per-mille with their result.
To get a first estimate of how low τ cut should be in order to reach asymptotic behavior for the τ cut → 0 extrapolation, we present in figure 3 the τ cut dependence of the NLO result computed with jettiness slicing. As before, the results are compared with the one obtained by using Catani-Seymour dipole subtraction that has already been implemented in MCFM. For the fitting to the asymptotic formula in eq. (2.1) we have used all data points. We have also limited the fitting data to sets of τ cut < 0.5 GeV and τ cut < 0.08 GeV. In both cases the fitted value of ∆σ NLO and its uncertainty stay the same, while of course increasing the fit uncertainty of the nuisance parameters affecting the tail. It is also possible to remove some of the smallest τ cut data points without altering the asymptotic value and uncertainty, to some extent. This makes us believe that the dominant behavior to extract ∆σ NLO is captured around the region τ cut 0.01 GeV.
Since our cross section uncertainties are not small enough to significantly fix the values of the nuisance parameters c 1 and Q in eq. (2.1), we do not include higher order power JHEP11 (2017)150 corrections. Including them deteriorates the fit quality of the nuisance parameters drastically. This follows the guidance in ref. [71]: in order to include additional terms in the fit, the model uncertainty induced by missing higher order power corrections in the fit should be comparable to the statistical uncertainties of the cross section values. For example the nuisance parameters in figure 3 are still fitted with standard uncertainties of 50-60%.
Adding the LO contribution to the NLO coefficient, our result of 71.09(1) fb compares well with the result from ref. [39] of 70.98 fb. They agree within 0.2%, but leave open the question of a possible tiny systematic difference in, for example, phase-space sampling. Our code for the NLO cross section with default settings is stable to values of τ cut at least as small as 10 −4 GeV, but at that point the numerical uncertainty using the same computational runtime increases drastically, of course. Instead we are able to obtain high precision results by sampling multiple larger values of τ cut using the additional information inferred from the asymptotic behavior in eq. (2.1). This also makes the result robust against statistical noise, underestimated uncertainties, or a possibly biased Monte Carlo grid, since our results for different values of τ cut are statistically independent. Note that even at τ cut = 0.02 GeV the deviation of the NLO coefficient from its asymptotic value is just about a percent. Using a value of τ cut 0.02 GeV for the calculation of the NNLO coefficient is then sufficient for any phenomenological application considering its relative contribution to the NLO cross section, as we will see.
We show the τ cut dependence of the NNLO cross section coefficient in figure 4. The NNLO coefficient extracted from the asymptotic fit contributes 8.3(1) fb to the full NNLO result of 80.7(1) fb. Our result is in excellent agreement with the result of 80.8(4) fb given in ref. [39]. We emphasize that the quoted uncertainty is purely statistical, originating from the Monte Carlo integration uncertainty in the weighted non-linear fit to the asymptotic τ cut behavior. Before continuing we first comment on the stability of our fit result, which we have rounded from 8.29(12) fb to 8.3(1). Removing the first and last point does not change the fitted asymptotic value within the uncertainty, as one obtains 8.25 (16) fb. The uncertainty induced by the actual fitting model is subleading, which is supported by the fact that adding for example a log 2 term does not change the fitted value and only increases the uncertainty slightly to 0.13 fb. Taking things to the extreme and using the NLO coefficient fitting form for our NNLO coefficient here, gives a result of 8.35(10) fb which again lies within our quoted uncertainty. The precision of our data points is not sufficient to resolve the additional logarithmic power and nuisance parameters, but instead relies on the general logarithmic shape.
For the charged lepton decay channel our LO result of 84.115(5) fb agrees with the result of 84.09 fb in table 3 of ref. [39] within 0.03%. Our NLO coefficient obtained from τ cut → 0 extrapolation is 53.89(4) fb and agrees perfectly with the C.-S. result of 53.92(1) fb within half a per-mille. This results in a total NLO cross section of 153. 15(2) fb that compares well with 153.1 fb given in table 3 of [39]. Using the fitted NNLO coefficient shown in figure 5 of 22.0(2) fb we obtain a total NNLO cross section of 178.2(2) fb, with a relative uncertainty of about one per-mille. We compare this to the value of (180 ± 1) fb in ref. [39] and find only broad compatibility. We note that their uncertainty is considerably larger than ours. This does not indicate a definite discrepancy but it is possible that the choice of slicing parameter could lead to such a difference. We observe that if we had instead chosen a larger value of τ cut ∼ 0.1 GeV then our results for both decay modes would have been compatible with those of ref. [39], within their quoted uncertainties. We also find similar behavior when comparing with the 7 TeV result. Our 7 TeV NLO calculation agrees within half per-mille with that reported in ref. [39]. The results of ref.
[39] predict a 7 TeV NNLO coefficient (using NNLO PDFs) of 19.6(9) fb. Our result for the same quantity is 18.3(2) fb and again only agrees marginally, suggesting a small systematic difference.

Standard Model phenomenology at 13 TeV
In order to present new phenomenological results and to study how the extrapolation of τ cut → 0 performs at a higher center of mass energy, we calculate cross sections and kinematical distributions for the neutrino decay channel at √ s = 13 TeV in the following section 4.1. Note again that we use different photon isolation cuts here compared to the 8 TeV comparison results in section 3; specifically we set = 0.1, n = 2 and R 0 = 0.4 (see paragraph 2). In section 4.2 we study the 13 TeV m llγ spectrum as a background to the very rare H → Z(→ ll)γ decay.

Decay to neutrinos
We expect that the cuts for the upcoming 13 TeV Zγ ATLAS analysis will be very similar to the ones that are summarized in table 3. Using these we investigate NNLO predictions for the cross section and kinematical distributions, where in the latter case we discuss  Table 3. Applied cuts for the Z →νν decay channel at a center of mass energy √ s = 13 TeV. the systematic uncertainty induced by using a finite value of τ cut instead of performing a systematic extrapolation procedure.
We must first re-establish the appropriate range of τ cut values that represents asymptotic behavior for our new set of cuts with √ s = 13 TeV. As before, we compute the τ cut dependence of the NLO cross section coefficient and compare it to the C.-S. obtained result. We find that the same range of τ cut leads to asymptotic behavior, and values of τ cut 0.1 GeV already lead to less than a percent difference from the extrapolated value. Again, our fitted value of 20.45(1) is in perfect agreement with the C.-S. result of 20.439(3) fb. We find a similar situation at NNLO, as illustrated in figure 6, which shows the τ cut dependence of the NNLO coefficient. Using the fitted extrapolated value of (2.24 ± 0.12) fb allows us to calculate the total NNLO cross section as (86.09 ± 0.12) fb.
Our results for the cross section under this set of cuts, at each order in QCD, are shown in [fb] Figure 6. τ cut dependence of the NNLO cross section coefficient in the neutrino decay channel at √ s = 13 TeV with cuts given in table 3. The solid black line is a fit to the expected asymptotic behavior in eq. (2.2). The solid red line shows the asymptotic value for τ cut → 0 of ∆σ NNLO = (2.24 ± 0.12) fb. The dashed red lines show the 95% uncertainty band of about 0.12 fb. tainty, obtained using the 7-point variation discussed in paragraph 2. We see that these scale uncertainties are reduced when going from NLO to NNLO but remain well above the residual 95% confidence level uncertainty that results from performing the fit to the τ cut → 0 limit. We therefore conclude that the use of the jettiness slicing method, and our extrapolation procedure, does not lead to any systematic effect of phenomenological relevance. In table 4 we also indicate the expectation for the jet-vetoed cross section (N jet = 0), although this prediction is susceptible to large logarithmic contributions that may require resummation, see for example a recent treatment of such issues for W + W − production [77].
We now turn to a study of NNLO predictions for kinematical distributions. We begin by considering a differential evaluation of the statistical and systematic numerical uncertainties given by integration precision and finite τ cut , just as we have already done for the total cross section. This is illustrated in figure 7, where we plot the ratio of the NNLO p γ For reference we present the scale uncertainty band obtained for τ cut = 0.018 GeV. The first observation from the plot is that the residual systematic uncertainty induced by a finite τ cut is small compared to the remaining integration uncertainty and the scale variation uncertainty. Additionally, the numerical integration uncertainty is small compared to the scale variation uncertainty. Although the scale uncertainty dominates the total, if it is imperative to obtain results with uncertainties of 0.1% then a more costly systematic extrapolation of τ cut → 0 can certainly be performed. This could be performed in a manner similar to the total cross section studies presented in table 4 and figure 6. For simplicity, since the systematic residual τ cut dependence is small compared to the scale uncertainty, henceforth all our studies of distributions will be performed using τ cut = 0.018 GeV.
In figures 8-10 we show results for the jet inclusive and 0-jet exclusive photon transverse momentum distributions, as well as the ννγ transverse mass distribution, respectively. Figure 8 shows the jet inclusive photon transverse momentum distribution at different perturbative orders as well as the k-factor relative to the NLO result. As anticipated in the discussion above, the numerical uncertainty of the NNLO results of about 0.2% is subleading with respect to the scale uncertainty band. The impact of the perturbative corrections for the new set of 13 TeV cuts is much smaller than for the 7 TeV ATLAS analysis. Whereas the corrections for the 7 TeV cuts are between 10 and 20% for p γ T < 250 GeV [39], the perturbative corrections for our set of cuts are between 4 −8% and are flat for large p γ T . In figure 9 we consider the photon transverse momentum distribution with a jet veto, where the NNLO corrections increase the cross section by up to about 10% at p γ T ∼ 1 TeV.  table 3. The upper panel shows the differential jet inclusive cross section, while the lower panel shows the ratio to the differential NLO distribution including a scale variation uncertainty band. The numerical integration uncertainty is 0.2% for the NNLO result. away from the production threshold. For these distributions we have also checked again that the systematic error due to a finite τ cut is small compared to the other uncertainties.

Decay to charged leptons as background to a Higgs signal
A particularly interesting aspect of the decay channel to electrons (charged leptons) is its importance as a background to Higgs production with subsequent decay to Zγ. In this subsection we consider this case using the cuts given in table 5, which are motivated by the 8 TeV ATLAS H → Zγ search in ref. [78]. Our cuts resemble those for the search in the e + e − decay channel. For the muon decay channel the cuts in ref. [78] are altered. ratio to NLO Figure 9. 13 TeV p γ T distribution for the neutrino decay channel with the cuts in table 3. The upper panel shows the differential jet vetoed cross section, while the lower panel shows the ratio to the differential NLO distribution including a scale variation uncertainty band. The numerical integration uncertainty is about 0.2 to 0.5%.
Since the H → Zγ signal in the charged lepton plus photon invariant mass spectrum m llγ will be tiny compared to the direct Zγ production background an understanding of the m llγ background shape will be essential. 6 This is despite the fact that experimental analyses generally rely on polynomial sideband fits. Nevertheless, perturbative corrections must be calculated to make sure that no unexpected shape corrections are introduced until the theory uncertainties are under good control at NNLO and beyond. To study this, we show the m llγ spectrum in figure 11 and focus on the NNLO to NLO k-factor: the JHEP11(2017)150 distribution for the neutrino decay channel decay with cuts given in table 3. The upper panel shows the differential cross section, while the lower panel shows the ratio to the differential NLO distribution including a scale variation uncertainty band. The numerical integration uncertainty is 0.2% for the NNLO results.
perturbative NNLO corrections are fortunately relatively flat and about 8-14% relative to NLO. Since the shape of the spectrum is stable as the perturbative order increases, it may be beneficial to compare the data (and any resulting polynomial fits) to a theoretical prediction such as this one. Deviations between the two could thus be used to constrain experimental issues such as mis-tagged background events.
We have also computed the theoretical expectation for the signal process at NNLO under the same set of cuts. Our predictions for the corresponding H → Zγ signal cross section in gluon fusion production in the infinite top-mass limit are given in table 6. As    is typical for gluon fusion Higgs production, the perturbative corrections are large and a big scale uncertainty remains at NNLO. To check whether the signal to background ratio changes noticeably when taking into account the NNLO corrections, we also include the background cross section in table 6, which has been obtained by adding the bins in the range 120 GeV to 130 GeV from figure 11. Depending on the experimental analysis this range might have to be increased, of course, further suppressing the S/B ratio. While the central S/B ratio of about four per-mille increases slightly towards NNLO from about 3.5 per-mille at NLO, the substantial scale uncertainties at NLO already account for such a change.

Anomalous couplings and probe for new physics
In this section we study anomalous ZZγ and Zγγ coupling contributions introduced by field operators up to dimension 8, requiring Lorentz invariance and electromagnetic gauge invariance [14,28,58,59]. The effective ZγZ vertex is described by where h Z i , i = 1, . . . , 4 are the effective ZγZ coupling factors and Λ is conventionally chosen to be the Z boson mass. For a different scale the coupling factors have to be scaled accordingly [14]. The couplings for i = 1, 2 are CP-violating, whereas for i = 3, 4 they preserve CP symmetry. The vertex for Zγγ can be obtained by setting q 2 1 = 0 and replacing h Z i by h γ i . Since the CP-conserving and CP-violating couplings do not interfere, their sensitivities to anomalous couplings are nearly identical and usually only h V 3 and h V 4 are considered in analyses.
The couplings h Z,γ i are zero at tree level in the SM and thus provide a unique opportunity to probe for new physics. We have implemented the vertices in our NNLO code JHEP11(2017)150 to allow for a unified calculation of NNLO Zγ production in the SM and in the presence of anomalous couplings. Whereas previously in experimental analyses anomalous coupling contributions were calculated at NLO using MCFM-8.0 and combined with a NNLO SM calculation [9], this error prone combination is no longer necessary with the additional benefit of having a consistent NNLO prediction.
In general the effect of anomalous couplings is to cause the matrix elements to grow at high energies so that it is this region that has the highest sensitivity to their presence. When comparing limits one has to be careful though, as limits obtained from different kinematic regions can be based on different assumptions on the scale of new physics Λ [79]. To partially cure this problem and to avoid unitarity violation, a form factor that dampens the anomalous coupling contribution at higher energies can be introduced. In our study we do not apply any such form factor.
Measurements of Zγ production at the LHC that focus on probing anomalous couplings have been performed by CMS and ATLAS at 7 and 8 TeV [1, 3-6, 8-10]. Previously anomalous couplings have been constrained by CDF and D0 at the Tevatron [80,81]. In their latest 8 TeV analysis the ATLAS collaboration derives limits both using a form factor with a scale of 4 TeV and no form factor [9], while the CMS collaboration uses no form factor [10]. In all cases the obtained limits are of the order 10 −3 on h V 3 and 10 −5 on h V 4 . These results set the scale for the appropriate ranges of anomalous couplings that we consider shortly.
For the current best limits from ATLAS and CMS, the ATLAS analyses use an exclusive zero-jet selection with cuts of p γ T > 250 GeV for the Z →ll decay, and p γ T > 400 GeV for the Z →νν decay channel. The neutrino channel has the highest sensitivity to anomalous couplings due to the larger branching fraction. ATLAS and CMS also give measured cross sections, which currently show a huge downward fluctuation with respect to the SM. Such a big negative contribution cannot be achieved or explained with the anomalous coupling contributions because the couplings need to be small enough to produce a negative interference with the SM that is not overwhelmed by the squared anomalous coupling term. Indeed, the cross section dependence on the anomalous couplings is almost purely quadratic, which explains the high symmetry of the experimental anomalous coupling limits about zero.
Since for the anomalous coupling extraction a jet veto and a high p T photon are required, we want to stress a caveat of QCD fixed order calculations in this region. For most precise predictions the jet veto requires resummation of its induced logarithmic terms, while the high p T photon mandates the inclusion of EW corrections [34][35][36]. Currently such jet veto logarithms have not been resummed for Zγ diboson production, but for example for W + W − production [77]. For color singlet Z production they have been studied [82] and can require a more conservative estimation of perturbative uncertainty if left out [83]. Our results show the impact of fixed order NNLO QCD corrections, and for a comparison with experimental data to extract limits on anomalous couplings it is important to take into account these additional resummed and EW contributions.
We do not want to exhaustively cover all anomalous couplings here since their behavior from a theoretical viewpoint is mostly similar. We thus give two representative examples and, in each case, consider ranges of anomalous couplings that directly cover the current limits discussed above. First we study the dependence of the cross section on just h Z 3 as a one-dimensional example of the impact of NNLO corrections on anomalous coupling limits. Figure 12 shows the effect of h Z 3 on the cross section in the neutrino decay channel at 13 TeV. We apply the 13 TeV neutrino channel cuts given in table 3 with an additional zero-jet veto and a photon transverse momentum requirement of p γ T > 400 GeV. The change in the central prediction due to the NNLO corrections is relatively small, around 5-10% over the whole range of h Z 3 , and negative. However, accounting for the effect of the entire scale uncertainty band allows for negative corrections of up to 20-40%. Even more interesting is to examine how these perturbative corrections translate into a change in the estimated limit on h Z 3 . To directly quantify the change in the estimated limit after accounting for the NNLO corrections we use the following procedure. We fit a quadratic function to describe the dependence of the cross section σ on the anomalous coupling h, σ(h) = σ SM + hσ interf. + JHEP11(2017)150  Figure 13. The change in the limit on the anomalous coupling h Z 3 when going from NLO to NNLO as explained in the text. The grey band is obtained by using the scale uncertainty in the theoretical prediction.
h 2 σ sq. . Here σ SM is the SM cross section and the terms proportional to h and h 2 represent the effects of interference with the anomalous coupling term and the anomalous coupling squared term, respectively. Inverting the equation gives the value of the anomalous coupling as a function of the cross section that is measured. We then plot the ratio of the NNLO prediction of σ NNLO (h) to the NLO prediction σ NLO (h). This number represents the amount the limit on an anomalous coupling is changed by the NNLO corrections. For example, for h Z 3 this amounts to a loosening of the limit by 10-15% using our the central scale choice, driven by the fact that the NNLO corrections are negative. We do the same with the envelope given by the scale variation in order to derive an uncertainty band for the effect on the limit. The result of this procedure is displayed in figure 13. Other anomalous couplings show a similar behavior, with the NNLO corrections leading to a similar impact on both the central value that can be extracted and the associated scale uncertainty.
Having outlined the procedure, we now turn to a two-dimensional analysis performed in the h Z 3 , h Z 4 plane. In figure 14 we show a contour plot of σ(h Z 3 , h Z 4 ) at LO, NLO and NNLO with contour lines assuming a cross section that is 1.1, 2 and 5 times the SM prediction at the corresponding perturbative order. The axis limits are chosen to roughly match the current 2D limits on h Z 3 , h Z 4 ; see figure 14 in ref. [9]. For the plot, we fit the calculated σ(h Z 3 , h Z 4 ) cross section at interpolation points to the functional form The 2D plot essentially confirms the findings of the one-dimensional analysis. The NNLO corrections to the cross section in the presence of anomalous couplings are relatively small (for the used central scale), as in the SM, but nevertheless they should be included in order to provide a consistent extraction of limits at NNLO.

Conclusions and outlook
We have presented a calculation of Zγ production at NNLO, fully including the neutrino and charged lepton decays, and in the latter case also including the virtual photon contribution. Since we are using the jettiness slicing method, the recent publication of power corrections [71,76] for Drell-Yan type color singlet processes motivated a study of their impact on our calculation. For Zγ production, in addition to the Drell-Yan s-channel type of contribution, also a t-channel diagram contributes, where the photon is not radiated from the final state charged leptons. Since the existing power corrections are calculated solely for the Drell-Yan type contribution their inclusion does not dramatically improve the τ cut dependence, unless specific cuts are employed. We implemented the boosted definition of the jettiness observable as the relevant slicing parameter and found the improvements to be moderate for the cuts used in Zγ analyses with relatively central Zγ systems. Since the improvement is consistent and comes essentially for free, we used this definition throughout our study.

JHEP11(2017)150
For the validation of our ingredients and the implementation itself, we compared with the NNLO results in ref. [39] and found good agreement in the neutrino decay channel but only limited agreement in the charged lepton decay channel. Our matrix elements have been checked with comparisons against other codes as described in the text. To obtain reliable estimates for our cross section predictions we used an extrapolation procedure of τ cut → 0 determined by the form of the leading jettiness power corrections. We note that if we had used a value of τ cut ∼ 0.1 GeV, we would have found agreement with the existing results for both channels, within uncertainties. The difference could be due to the value of the slicing parameter in the previous calculation.
For 13 TeV SM phenomenology, we studied Zγ in the neutrino decay channel in differential distributions with cuts similar to an upcoming ATLAS analysis. We have shown that by using a sufficiently low value of τ cut , the systematic error made is small compared to the numerical integration uncertainty in phenomenological applications, which again is small compared to the residual scale uncertainty. For the charged lepton decay channel we envisioned our direct Zγ production as a background to H → Zγ production. We studied the invariant mass distribution of the Zγ system and have shown that NNLO effects are relatively flat at about 8-14% between the Higgs mass and 165 GeV.
Finally, we presented results of our implementation of the Zγγ and ZZγ anomalous couplings at NNLO. We considered a one-dimensional example and a two-dimensional contour plot in the region of current experimental limits. The perturbative corrections for the jet-vetoed set of cuts with p γ T > 400 GeV are about 5-10% and will allow for a more precise extraction of anomalous coupling limits. Our implementation, which will be released publicly in MCFM, will allow for a unified NNLO SM and anomalous couplings analysis.
Diboson production, and specifically Zγ production, belongs to a class of processes with large perturbative corrections and severely underestimated scale uncertainties at NLO. This is due to the fact that up to NNLO new partonic channels open up. One thus hopes that for the first time with the inclusion of NNLO corrections the scale uncertainty estimate is reliable, even though the new channels at NNLO are only in leading order in their partonic sense. Including higher order effects, for example for the gg channel as has been done in diphoton production [41], is possible since all ingredients are available. A full N 3 LO calculation could give a definite answer and should be on theorists' agenda anyway, since the uncertainties of results from ATLAS and CMS are already competing with the estimated NNLO theory precision.

JHEP11(2017)150
publisher, by accepting the article for publication, acknowledges that the U.S. Government retains a non-exclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this manuscript, or allow others to do so, for U.S. Government purposes.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.