W γ and Zγ production at the LHC in NNLO QCD

: We consider the production of Wγ and Zγ pairs at hadron colliders. We report on the complete fully differential computation of radiative corrections at next-to-next-to-leading order (NNLO) in QCD perturbation theory. The calculation includes the leptonic decay of the vector boson with the corresponding spin correlations, off shell effects and final-state photon radiation. We present numerical results for pp collisions at 7 and 8 TeV and we compare them with available LHC data. In the case of Zγ production, the impact of NNLO corrections is generally moderate, ranging from 8% to 18%, depending on the applied cuts. In the case of Wγ production, the NNLO effects are more important, and range from 19% to 26%, thereby improving the agreement of the theoretical predictions with the data. As expected, the impact of QCD radiative corrections is significantly reduced when a jet veto is applied. Abstract: We consider the production of W γ and Zγ pairs at hadron colliders. We report on the complete fully diﬀerential computation of radiative corrections at next-to-next-to-leading order (NNLO) in QCD perturbation theory. The calculation includes the leptonic decay of the vector boson with the corresponding spin correlations, oﬀ shell eﬀects and ﬁnal-state photon radiation. We present numerical results for pp collisions at 7 and 8 TeV and we compare them with available LHC data. In the case of Zγ production, the impact of NNLO corrections is generally moderate, ranging from 8% to 18%, depending on the applied cuts. In the case of W γ production, the NNLO eﬀects are more important, and range from 19% to 26%, thereby improving the agreement of the theoretical predictions with the data. As expected, the impact of QCD radiative corrections is signiﬁcantly reduced when a jet veto is applied. on


Introduction
The discovery of a new scalar resonance in the search for the Standard Model (SM) Higgs boson [1,2] is a milestone in the LHC physics programme. The properties of this new particle closely resemble those of the Higgs boson, but further work is needed to clarify if it is really the Higgs boson predicted by the SM, or something (slightly) different. Vectorboson pair production has a prominent role in this context. It represents an irreducible background to Higgs and new-physics searches, and, at the same time, it provides information on the form and the strength of the vector-boson gauge couplings. The interactions of W and Z bosons with photons are particularly interesting as they test the W W γ and ZZγ couplings, which are predicted by the non-Abelian SU(2) L ⊗ U(1) Y gauge group. Constraints on W W γ and ZZγ anomalous couplings have been obtained at LEP [3]. At hadron colliders, studies of V γ final states have been first carried out at the Tevatron [4][5][6], and they were used to set more stringent limits on anomalous couplings. The high-energy proton-proton collisions at the LHC allow us to explore the production of V γ (V = W ± , Z) pairs in a new energy domain. Measurements of V γ final states have been carried out by ATLAS [7][8][9][10] and CMS [11][12][13][14] using the data sets at centre-of-mass energy √ s = 7 and 8 TeV. These measurements have been compared to the SM predictions and used to improve the limits on anomalous couplings and on the production of possible new resonances.

JHEP07(2015)085
When considering the V γ final state, besides the direct production in the hard subprocess, the photon can also be produced through the fragmentation of a QCD parton, and the evaluation of the ensuing contribution to the cross section requires the knowledge of a non-perturbative photon fragmentation function, which typically has large uncertainties. The fragmentation contribution is significantly suppressed by the photon isolation criteria that are necessarily applied in hadron-collider experiments in order to suppress the large backgrounds. The standard cone isolation, which is the standard choice in the experiments, suppresses a large fraction of the fragmentation component. The smooth cone isolation completely suppresses the fragmentation contribution [15], but the algorithm is difficult to be implemented experimentally.
The present status of theoretical predictions for V γ production at hadron colliders is as follows. The V γ cross section is known in next-lo-leading-order (NLO) QCD [16,17], and the leptonic decay of the vector boson has been included in ref. [18]. In the case of Zγ the loop-induced gluon fusion contribution, which is formally next-to-next-to-leading order (NNLO), has been computed in refs. [19,20], and the leptonic decay of the Z boson, together with the gluon-induced tree-level NNLO contributions from gg → Zγ qq, have been added in ref. [21]. The NLO calculation for V γ, including photon radiation from the final-state leptons, the loop-induced gluon contribution and the photon fragmentation at LO have been implemented into the general purpose numerical program MCFM [22]. Electroweak (EW) corrections to V γ production have been computed in refs. [23,24]. The full NLO EW corrections to W γ production with leptonically decaying W bosons, taking into account all off-shell effects in the complex-mass scheme, and all effects originating from initial-state photons, have been computed in ref. [25]. For W γ production, the NLO computation has been matched to a parton shower according to the MiNLO prescription [26] in ref. [27].
In a previous letter [28] we have presented the results of the first complete NNLO calculation for Zγ production. In the present paper, we extend this calculation to the complete class of V γ production with leptonic decays, namely to both Zγ production with visible (Z → ℓ + ℓ − ) and invisible (Z → νν) Z-boson decays, and to W γ production with the respective decays W + → νℓ + and W − → ℓ −ν . Off-shell effects and final-state photon radiation are consistently included. 1 For these production channels, we present detailed results on fiducial cross sections and distributions at √ s = 7 and 8 TeV, and provide comparisons to ATLAS data, where available. The paper is organized as follows. In section 2 we provide some technical details of our computation and discuss the particular challenges in the cancellation of infrared singularities. Section 3 contains our theoretical predictions for all V γ processes as well as a comparison with experimental data. We consider Zγ production in the visible (Z → ℓ + ℓ − ) and invisible (Z → νν) decay channels in section 3.3 and section 3.4 and W γ production in section 3.5. In section 3.6 we discuss the different impact of QCD radiative corrections in the W γ and Zγ processes and its physical origin. In section 4 we summarize our results and comment on the remaining uncertainties.

Details of the calculation
In this section we discuss the details of our calculation. We first point out that the notation "V γ" suggests the production of an on-shell vector boson plus a photon, followed by a factorized decay of the vector boson. Instead, we actually compute the NNLO corrections to the processes pp → ℓ + ℓ − γ +X, pp → ν ℓνℓ γ +X, and pp → ℓν ℓ γ +X, where, in the first case, the lepton pair ℓ + ℓ − is produced either by a Z boson or a virtual photon. All contributions where the final-state photon is radiated off¡ the charged leptons are consistently included (see figure 1 and figure 2). The shortcuts "Zγ" and "W γ" are used only for convenience.
The NNLO computation requires the evaluation of tree-level scattering amplitudes with up to two additional (unresolved) partons, of one-loop amplitudes with up to one additional (unresolved) parton [30,31], and of one-loop squared and two-loop corrections to the Born subprocess (qq → ℓ + ℓ − γ and qq → ν ℓνℓ γ for Zγ, qq ′ → ℓν ℓ γ for W γ). Furthermore, processes with charge-neutral final states receive loop-induced contributions from the gluon fusion channel (gg → ℓ + ℓ − γ and gg → ν ℓνℓ γ). In our computation, all required treelevel and one-loop amplitudes are obtained from the OpenLoops generator [32], 2 which implements a fast numerical recursion for the calculation of NLO scattering amplitudes within the SM. For the numerically stable evaluation of tensor integrals we rely on the Collier library [33], which is based on the Denner-Dittmaier reduction techniques [34,35] and the scalar integrals of [36].
The two-loop corrections to the Drell-Yan-like Born processes, where the photon is radiated off the final-state leptons, have been available for a long time [37]. The last missing ingredient, the genuine two-loop corrections to the V γ amplitudes, have been presented in ref. [38].
In section 2.1, we give a sketch of the q T subtraction formalism, the method we use to combine the (separately divergent) ingredients of the NNLO calculation to obtain quan-JHEP07(2015)085 titative predictions. The bookkeeping of all partonic subprocesses and the numerical integration of the different cross section contributions is managed by the fully automatized Munich framework, 3 which is described in section 2.2. In section 2.3, the double-virtual and gluon fusion contributions are discussed, and section 2.4 addresses the real-emission part of the NNLO cross section with a discussion of the numerical stability issues, namely the reliable evaluation of one-loop amplitudes in the real-virtual part and the numerically stable phase-space integration in the double-real part.

The q T subtraction formalism
The implementation of the various scattering amplitudes in a complete NNLO calculation is a highly non-trivial task due to the presence of infrared (IR) singularities at intermediate stages of the calculation, which prevents a straightforward implementation of numerical techniques. Various methods have been proposed and used to overcome this difficulty. They are based either on extensions of the subtraction method [39][40][41][42] at NNLO [43][44][45][46], or on sector decomposition [47,48]. More recently, also a combination of the subtraction method with sector decomposition has been proposed [49,50]. The q T subtraction formalism [51] is an independent method to handle and cancel IR singularities at the NNLO. In its present formulation the method applies to the production of a colourless high-mass system F in generic hadron collisions, and has been applied to the computation of NNLO corrections to several hadronic processes [28,[51][52][53][54][55][56][57]. According to the q T subtraction method [51], the pp → F + X cross section at (N)NLO can be written as where dσ F+jet (N)LO represents the cross section for the production of the system F plus one jet at (N)LO accuracy, and can be evaluated with any available NLO subtraction formalism. The (IR subtraction) counterterm dσ CT (N)NLO is obtained from the resummation of logarithmically-enhanced contributions to q T distributions [58]. The practical implementation of the contributions in the square bracket in eq. (2.1) is described in more detail in section 2.4.
The 'coefficient' H F (N)NLO encodes the loop corrections to the Born-level process and also compensates 4 for the subtraction of dσ CT (N)NLO . It is obtained from the (N)NLO truncation of the process-dependent perturbative function The NLO calculation of dσ F requires the knowledge of H F(1) , and the NNLO calculation also requires H F (2) . The general structure of H F(1) is known [59]: H F(1) is obtained 3 Munich is the abbreviation of "MUlti-chaNnel Integrator at Swiss (CH) precision" -an automated parton level NLO generator by S. Kallweit. In preparation. 4 More precisely, while the behavior of dσ CT (N)NLO as qT → 0 is dictated by the singular structure of dσ F+jet (N)LO , its non-divergent part in the same limit is to some extent arbitrary, and its choice determines the explicit form of H F (N)NLO .

JHEP07(2015)085
from the process-dependent scattering amplitudes by using a process-independent relation. Exploiting the explicit results of H F (2) for Higgs [60] and vector boson [61] production, the process-independent relation of ref. [59] has been extended to the calculation of the NNLO coefficient H F(2) [62]. These results have been confirmed through an independent calculation in the framework of Soft-Collinear Effective Theory (SCET) [63,64]. The counterterm dσ CT (N)NLO only depends on H F (N)LO , i.e. for a NNLO computation, it requires only H F(1) as input, which can be derived from the one-loop amplitudes to the involved Born subprocesses.
As can be seen from the discussion above, a significant part of the ingredients needed to perform the NNLO computation are actually NLO-like in nature, allowing us to build the implementation of our NNLO calculation upon existing NLO tools. We have based this calculation on the Munich framework, which is briefly discussed in the following. Necessary extensions of this framework are addressed as well in the rest of this section.

Organization of the calculation within the Munich framework
Munich is a fully automatized framework for the computation of fixed-order cross sections for arbitrary SM processes up to NLO accuracy, written in C++. After the hadronic process has been specified, Munich takes care of the bookkeeping of all partonic subprocesses that need to be taken into account. It automatically generates adequate phase-space parametrizations for each partonic subprocess by exploiting the resonance structure of the underlying (squared) tree-level Feynman diagrams. These parametrizations are combined using a multi-channel approach to simultaneously flatten the, in general complicated, resonance structure of the amplitudes and thus guarantee a reasonable convergence of the numerical integration. Several improvements like an adaptive weight-optimization procedure are implemented as well.
Munich was originally developed for the NLO calculations of [65][66][67], where only massless colour-charged particles were involved. To account for the mediation of IR singularities between the different phase spaces of virtual and real corrections, the Catani-Seymour dipole subtraction method [41,42] was implemented. In ref. [68], the framework was extended to massive QCD according to ref. [69], and in ref. [70] the extension to complete SM calculations, i.e. including also EW corrections at NLO, was presented. All dipole terms necessary to numerically cancel the soft and collinear divergences of the real corrections are constructed at runtime from spin-and colour-correlated matrix elements of the underlying partonic Born subprocesses. Moreover, additional phase-space parametrizations based on the reduced dipole kinematics are generated and included in the multi-channel, supplementing the parametrizations based on real-emission kinematics. Analogously, the analytically integrated dipoles, the so-called K +P terms and the I-operator, which compensate for the dipole terms subtracted on the real-emission side, are automatically constructed at runtime from colour-correlated matrix elements of the Born subprocesses. For the evaluation of all involved matrix elements up to the one-loop order, Munich provides an automatic interface to amplitudes generated by the OpenLoops generator [32].
The guiding principle in our NNLO implementation is to keep the additionally introduced process dependence to the bare minimum, essentially limiting it to the two-loop JHEP07(2015)085 amplitudes entering H F in eq. (2.1). All other process dependent information entering the various pieces in eq. (2.1) have been expressed in terms of NLO quantities available inside Munich via its interface to OpenLoops.

Double-virtual and gluon fusion contribution
In our implementation of the q T subtraction formalism, the only component that needs to be provided on a process-by-process basis are the two-loop -and one-loop-squared 5amplitudes to the Born processes, which enter the coefficient H F (N)NLO in eq. (2.1). Ref. [38] provides the analytical expressions of the resonant helicity amplitudes for V γ production (e.g. qq → Zγ → ℓ + ℓ − γ) in terms of two-dimensional harmonic polylogarithms up to weight 4. The non-resonant (final-state radiation) contribution (e.g. qq → ℓ + ℓ − → ℓ + ℓ − γ) is described by the two-loop quark form factor from ref. [37]. We have implemented the helicity amplitudes directly into a C++ library. For the numerical evaluation of the harmonic polylogarithms we use the tdhpl library [71]. Our implementation allows for several thousand amplitude evaluations per minute, rendering the time spent on computing the two-loop contribution negligible compared to the time needed for the real emission contribution.
For the processes with an electrically neutral final state at Born level, i.e. pp → ℓ + ℓ − γ and pp → ν ℓ ν ℓ γ, the loop-induced gluon fusion processes gg → ℓ + ℓ − γ and gg → ν ℓ ν ℓ γ represent a separately IR-finite and gauge-invariant part of the cross section. Though of O α 2 S , and thus formally of NNLO, the gluon fusion contribution is often included already at NLO (see for example MCFM [22].) The reason for that is the assumption that the large gluon luminosities at the LHC could compensate for the additional α S suppression, and hence the contribution could numerically become as important as the NLO corrections themselves. As we will see in section 3.3 and section 3.4, this is not the case in Zγ production. The gluon fusion contribution is small even when compared to the other NNLO corrections. In our computation, these (finite) one-loop-squared gluon fusion amplitudes are evaluated using OpenLoops.

Real-emission and counterterm contribution
All NNLO contributions with vanishing total transverse momentum q T of the final state system F are collected in the coefficient H F NNLO , which is discussed in the previous section. The remaining part of the NNLO cross section, namely the difference in the square bracket in eq. (2.1), is formally finite in the limit q T → 0, but the terms separately exhibit logarithmic divergences in this limit. The counterterm dσ CT (N)NLO is integrated over the n-particle Born phase-space, 6 while the term dσ F+jet NLO in eq. (2.1), called the real-emission contribution in the following, involves two different phase spaces with one and two additional QCD partons, respectively, and thus needs to be integrated separately over the respective (n + 1)and (n + 2)-particle phase spaces. To achieve a numerical cancellation of the singularity, a technical cut on q F T needs to be introduced to render both terms separately finite. In this 5 The OpenLoops generator also provides the non-finite one-loop-squared matrix elements in terms of the corresponding coefficients of the Laurent series. These results are, however, only used as a numerical check in this calculation. 6 The integration over qT only acts on an explicit qT dependence of the phase-space weight.

JHEP07(2015)085
sense, the q T subtraction method works very similarly to a phase-space slicing method at NLO. In practice, a technical cut, r cut , on the dimensionless quantity r ≡ q T /M , where M denotes the invariant mass of the final-state system F , turns out to be a more convenient choice than a cut on q T itself. By construction, the counterterm dσ CT (N)NLO cancels all q T -divergent (logarithmic) terms from the real-emission contributions, implying that the r cut dependence of their difference does not only vanish in the limit r cut → 0, but should become arbitrarily small for sufficiently small values of r cut . In practice, however, as both the counterterm and the realemission contribution grow arbitrarily large for r cut → 0, the statistical accuracy of the Monte Carlo integration degrades, preventing one from pushing r cut arbitrarily low. In general, the absence of any strong residual r cut dependence in the difference between the real contribution dσ V γ+jet (N)LO and the counterterm dσ CT (N)NLO provides a strong check on the correctness of the computation since any (significant) mismatch between the contributions would result in a divergent cross section in the limit r cut → 0. To monitor the r cut dependence without the need of several CPU-intensive runs, our implementation allows for simultaneous cross section evaluations at arbitrary r cut values.
Typical values of r cut are of the order of 10 −3 , requiring a numerically stable evaluation of dσ V γ+jet NLO down to the IR-divergent region where the transverse-momentum of the V γ system reaches values of about 0.1 GeV. While the IR-enhanced phase-phase region with a very-low-q T jet is excluded in standard NLO calculations, it needs to be resolved at very high precision in a NNLO computation. This issue obviously challenges a dedicated NLO program like Munich. We have slightly modified the phase-space generation to achieve the required precision and reliability of our NNLO results within the Munich framework.
The real-emission contribution to the NNLO cross section is, apart from the presence of a very-low-q T jet discussed above, a usual NLO calculation, and can thus be treated with the Catani-Seymour dipole subtraction method implemented within Munich. Consequently, it can be split into a real-virtual (RV), a real-collinear (RC) and a double-real (RR) contribution, Here, the mediation of all IR divergences (given a finite r cut value) by dipole subtraction is implicitly understood, i.e. each contribution on the right-hand side is finite and can be numerically integrated over the respective phase space. Whereas the real-collinear subcontribution does not exhibit any peculiar issues, the other two subcontributions on the right-hand side pose different challenges when integrated into the deep IR (q T 0) region. The real-virtual contribution requires the evaluation of one-loop matrix elements to V γ +jet production, which has become a standard task for automatic one-loop tools. However, as we evaluate the matrix elements in the deep IR region, i.e. far away from the phase-space region that is relevant for 2 → 3 hard scattering processes at NLO, numerical instabilities in the amplitude evaluation might be a source of concern. To address this issue, OpenLoops implements a fully automated system that monitors the numerical accuracy of loop amplitudes and cures possible instabilities at runtime. This system exploits the fact that, for the reduction to scalar integrals, OpenLoops allows one to flexibly switch from tensor-reduction [34,35] to OPP reduction [72] algorithms.

JHEP07(2015)085
To perform the reduction from tensor to scalar integrals, we make use of the Collier library [33] that consists of two independent implementations of the Denner-Dittmaier reduction algorithm [34,35]. The presence of potential residual instabilities is tested by comparing the two implementations of tensor reduction. In presence of instabilities that exceed one permille of the Born amplitude, the one-loop amplitude is reevaluated using the CutTools [73] implementation of OPP reduction in quadruple precision. In this case OneLOop [74] is used for the evaluation of scalar integrals. Finally, the accuracy of the result in quadruple precision is verified by a consistency check based on the rescaling of all dimensionful variables. For all processes considered in this paper we find that the Denner-Dittmaier reduction works fast and reliably for more than 99% of the phase space points. This allows one to restrict the usage of quadruple precision to a tiny fraction of points with a minor impact on the total runtime of the calculation.
The double-real contribution, on the other hand, only involves tree-level amplitudes, but contains an additional unresolved parton in the final state. Moreover, it involves several dipole terms located on different (n + 1)-particle phace spaces. Like any other phase-space cut, the restriction r > r cut has to be applied on the respective phase space and can thus lead to miscancellation issues, which are well-known from NLO calculations, but might be amplified here when happening in the deep IR region. These complications render the numerical phase-space integration in the deep IR region more challenging compared to the integration of the single-emission phase space. Munich already provides a state-of-the-art multi-channel phase-space integrator, which greatly speeds up computations compared to a classical VEGAS integration procedure, and guarantees a reasonably stable and reliable convergence behaviour. In the deep IR regions, however, even an advanced multi-channel approach might fail to capture all relevant contributions. We have thus extended the Munich integrator by an additional VEGAS-like importance sampling on top of the multichannel parametrization. This hybrid approach results in an extremely efficient and reliable phase-space integration.
The numerical information on the r cut dependence of the cross section is used to perform an extrapolation to r cut = 0, which can in turn be used to quantify the uncertainty due to a finite r cut value. This estimated uncertainty is combined with the usual statistical integration error to provide an overall estimate of the numerical precision of our NNLO prediction. For the processes considered in this paper, the extrapolation to r cut = 0 turns out to be non-trivial, due to the interplay with the photon isolation, which requires r cut to be pushed much lower than in other applications of q T subtraction [51][52][53][55][56][57]. In this respect, the most problematic process in this paper is W γ production, where the qg andqg channels, which are the only channels exhibiting the quark-photon collinear singularity at NLO, receive significant corrections at the NNLO. Figure 3 shows the r cut dependence for W + γ production at 7 TeV, both at NLO and at NNLO. The comparison to the r cut -independent NLO result based on Catani-Seymour subtraction indicates convincing agreement in the limit r cut → 0, but unveals per-cent level deviations for r cut ≈ 1%, thereby enforcing us to go to much lower r cut values also at NNLO. 7 Nonetheless, the procedure allows us to control integrated NNLO cross sections to better than 1%.

JHEP07(2015)085
pp → e + νeγ + X @ 7 TeV pT,γ > 15 GeV pp → e + νeγ + X @ 7 TeV pT,γ > 15 GeV Figure 3. r cut dependence of the cross-section predictions to a sample setup for W + γ production at NLO (left) and NNLO (right), together with the 1σ integration error. (Note that the results for different r cut values are not statistically independent.) The left plot additionally shows the (r cutindependent) Catani-Seymour subtraction based NLO result (red). The right plot also depicts the estimated error band (green) around the extrapolated result.

Setup
For the electroweak couplings we use the so-called G µ scheme, where the input parameters are G F , m W , m Z . In particular we use the values G F = 1.16639 × 10 −5 GeV −2 , m W = 80.399 GeV, m Z = 91.1876 GeV, Γ Z = 2.4952 GeV and Γ W = 2.1054 GeV. We set the CKM matrix to unity. We use the MMHT 2014 [75] sets of parton distribution functions (PDFs), with densities and α S evaluated at each corresponding order (i.e., we use (n + 1)loop α S at N n LO, with n = 0, 1, 2), and we consider N f = 5 massless quarks/antiquarks and gluons in the initial state.
The default renormalization (µ R ) and factorization (µ F ) scales are set to An estimate of missing higher-order contributions is obtained by performing µ F and µ R variations by a factor of two around the central value. We find substantial cancellations between renormalization and factorization scale variations in some of the calculations we are going to present if the restriction µ R = µ F is imposed. These cancellations are assumed to be purely accidental. To accomodate this well-known feature, we consider also antipodal variations of the two scales [22], i.e. setting µ R = ξµ 0 , µ F = µ 0 /ξ and varying ξ between 1 2 and 2. In summary, we estimate scale uncertainties by varying µ F and µ R simultaneously and independently in the range 0.5µ 0 and 2µ 0 with no constraint on the ratio µ F /µ R .
The present formulation of the q T subtraction formalism [51] is limited to the production of colourless systems F and, hence, it does not allow us to deal with the parton fragmentation subprocesses. Therefore, we consider only direct photons, and we rely on the smooth cone isolation criterion [15]. Considering a cone of radius r = (∆η) 2 + (∆φ) 2 around the photon, we require that the total amount of hadronic (partonic) transverse JHEP07(2015)085 energy E T inside the cone is smaller than E max T (r), where p γ T is the photon transverse momentum; the isolation criterion E T < E max T (r) has to be fulfilled for all cones with r ≤ R. All results presented in this paper are obtained with ǫ γ = 0.5, n = 1 and R = 0.4.

Comparison to experimental data
The smooth cone isolation prescription adopted in our calculation is not yet implemented in experimental analyses. Measurements are usually performed by using a fixed cone isolation prescription (which corresponds to eq. (3.1) with n = 0), and thus our isolation prescription does not exactly correspond to what is done in the experiment. However, the parameters ǫ γ and R needed to specify the smooth cone have natural counterparts in the definition of the fixed cone, while the precise choice of the smoothing function (in our case parametrized by n) does only have a mild impact on the final result. Furthermore, recent studies carried out in diphoton production [76] show that for sufficiently tight isolation parameters, smooth and hard cone isolation yield very similar results. For the processes in this paper, we verified at NLO that the difference between using smooth and hard cone isolation is at the 1 − 2% level, 8 i.e. well below the current experimental uncertainties and still smaller than the remaining theoretical uncertainties. We can thus safely compare our theoretical predictions with experimental data.
Since the first results of our work have appeared [28,29], we have provided numerical predictions for Zγ production to the CMS collaboration [14], and for W γ production to the ATLAS collaboration [77], as a background in the H → W W analysis. These predictions were obtained by using the experimental cuts adopted in the corresponding analyses. In the present paper, we limit ourselves to compare our predictions to the ATLAS results for W γ and Zγ at 7 TeV [9].

pp
In our calculation of pp → ℓ + ℓ − γ at √ s = 7 and 8 TeV we adopt the selection cuts used by the ATLAS collaboration [9], summarized in table 1. We require the photon to have a transverse momentum of p γ T > 15 GeV (soft p γ T cut) or p γ T > 40 GeV (hard p γ T cut) and pseudorapidity |η γ | < 2.37. Each of the charged leptons is required to have p ℓ T > 25 GeV and |η ℓ | < 2.47, and the invariant mass of the lepton pair must fulfill m ℓ + ℓ − > 40 GeV. We require the separation in rapidity and azimuth ∆R between the leptons and the photon to be ∆R(ℓ, γ) > 0.7. Jets are reconstructed with the anti-k T algorithm [78] with radius parameter D = 0.4. A jet must have p jet T > 30 GeV and |η jet | < 4.4. We require the separation ∆R between the leptons (photon) and the jets to be ∆R(ℓ/γ, jet) > 0.3. At √ s = 8 TeV, the jet definition is slightly adjusted by using |η jet | < 4.5 instead of |η jet | < 4.4,

JHEP07(2015)085
Frixione isolation with ε γ = 0.5, R = 0.4, n = 1 Jets anti-k T algorithm with D = 0.4 adapting to the ATLAS Run II standard. With respect to resolved jets in the final states, we will consider both the inclusive (N jet ≥ 0) and the exclusive (N jet = 0) case.
The predicted cross sections with the soft p γ T cut, including the theoretical uncertainties from scale variations obtained as described at the beginning of section 3.1, can be found in table 2. When going from NLO to NNLO the cross section increases by 8% (3%) in the inclusive (exclusive) case, respectively. The fiducial cross sections measured by ATLAS at 7 TeV [9] are also reported in table 2. Both the NLO and NNLO predictions are in agreement with the experimental result, and the NNLO corrections improve the agreement, especially in the inclusive case.
The reduced impact of QCD radiative corrections when going from the inclusive (N jet ≥ 0) to the exclusive (N jet = 0) case is a well known feature in perturbative QCD calculations [79]. A stringent veto on the radiation recoiling against the Zγ system tends to unbalance the cancellation between positive real and negative virtual contributions, possibly leading to large logarithmic terms. The resummation of these logarithmic contributions has been the subject of intense theoretical studies [80][81][82], especially in the important case of Higgs boson production. The reduced impact of radiative effects in the presence of a jet veto is often accompanied by a reduction of scale uncertainties. In the present case, since we are considering a process initiated by quark-antiquark scattering, the impact of radiative corrections is smaller than in Higgs boson production. However, a reduction of scale uncertainties from the N jet ≥ 0 to the N jet = 0 case is already visible in table 2, and may signal the need of more sophisticated (conservative) methods to estimate perturbative uncertainties [80,83].  Table 2. pp → ℓ + ℓ − γ cross sections with the soft p γ T cut (p γ T > 15 GeV). Scale uncertainties are obtained from independent variations of µ R and µ F around the central scale µ 0 , as described in section 3.1. The numerical uncertainty of the NNLO prediction from statistical error and finite r cut are conservatively estimated to be about 0.3%. The last column provides the corresponding results by ATLAS.  Table 3. pp → ℓ + ℓ − γ cross sections with the hard p γ T cut (p γ T > 40 GeV). Scale uncertainties are computed as in table 2. The numerical uncertainty of the NNLO prediction from statistical error and finite r cut is conservatively estimated to be about 0.6%.

JHEP07(2015)085
Beyond the cross section in the fiducial region, ATLAS has also provided the measured cross section differential in the photon transverse momentum. A comparison of the resulting distribution with our theoretical NLO and NNLO predictions is displayed in figure 4 for both the inclusive and the exclusive case. In particular at transverse momenta p γ T ≤ 100 GeV, the inclusion of NNLO corrections tends to improve the agreement between data and theory. The comparison of the theoretical predictions to the data in the high transversemomentum region p γ T > 100 GeV is quite delicate. First, the experimental uncertainty in this region is quite large. Then, EW corrections are expected to become sizable and negative due to large Sudakov logarithmic contributions [23,24].
In figure 5 we compare the NLO and NNLO predictions for the invariant-mass distribution of the ℓ + ℓ − γ system with the distribution provided by ATLAS in ref. [9]. For this measurement, ATLAS increases the transverse momentum cut on the photon to p γ T > 40 GeV: the corresponding cross sections are reported in table 3. The relative impact of radiative corrections is 79% and 17% when going from LO to NLO and from NLO to NNLO, respectively. We conclude that the corrections are significantly larger compared to the case in which the soft p γ T cut (p γ T > 15 GeV) is applied. As the m ℓ + ℓ − γ differential cross section in figure 5 is normalized by the fiducial cross section, sizeable NNLO corrections are visible only in the first bin, where the agreement with data is slightly improved. This implies that the NNLO/NLO ratio is almost constant for larger invariant masses.   The more pronounced higher-order corrections in the case in which a hard p γ T cut is applied can be understood by studying the ℓ + ℓ − γ invariant mass distribution in a finer binning, which is shown in figure 6 for both the soft and the hard p γ T cuts. When the soft p γ T cut is applied, the relative impact of the NNLO corrections is small in the region around the Z → ℓ + ℓ − γ peak, where the fiducial cross section receives its dominant contribution, and then slowly increases with the invariant mass. When the hard p γ T cut is applied, the Z → ℓ + ℓ − γ peak is not populated at all at LO as the applied cuts produce a lower bound at m ℓ + ℓ − γ ≈ 97 GeV in LO kinematics. The region below the boundary contributes sizably to the cross section, but is only populated beyond LO, i.e. in this region the NLO computation provides actually only a LO prediction. Hence the NNLO predictions effectively correspond  Figure 6. Invariant mass distribution of the ℓ + ℓ − γ system at LO (blue, dotted), NLO (red, dashed) and NNLO (green, solid) for the setup with p γ T > 15 GeV (left) and the setup with p γ T > 40 GeV (right). The loop-induced gluon fusion contribution is also shown (pink, dash-dotted). The lower panel shows the NNLO/NLO ratio, and the bands indicate theoretical uncertainty estimates from scale variations.
to the first perturbative correction, with a comparably large K factor of about 1.4. The lower bound on m ℓ + ℓ − γ for LO kinematics also exists with the soft p γ T cut, namely at m ℓ + ℓ − γ ≈ 66 GeV. However, in this case the Z → ℓ + ℓ − γ peak is populated already at LO, and the region below the cut does not significantly affect the fiducial cross section. As already expected from figure 5, the NNLO/NLO ratio in the hard p γ T case is almost independent of m ℓ + ℓ − γ above m ℓ + ℓ − γ ≈ 140 GeV. Figure 6 also shows the contribution from the loop-induced gluon fusion process, which, as explained in section 2, respresents a finite and gauge invariant subcontribution to the full NNLO result. This contribution is often argued to be potentially sizable due to the large gluon luminosities at the LHC. In our calculation, however, the gluon fusion contribution turns out to be small. It amounts only to around 6(9)% of the full O α 2 S correction and, correspondingly, to less than 1(2)% of the total fiducial cross section in the soft and the hard p γ T case, respectively.

pp → ν ℓ ν ℓ γ
In the pp → Zγ → ν ℓ ν ℓ γ analysis for proton-proton collisions at √ s = 7 TeV, we again use the selection criteria applied by ATLAS [9]: compared to the pp → ℓ + ℓ − γ analysis, the transverse-momentum cut on the photon is made harder (p γ T > 100 GeV), and a cut on the missing transverse momentum, i.e. the vectorial sum of the neutrino momenta, p νν T > 90 GeV, is imposed. The jet algorithm, the photon isolation and all other eventselection criteria are the same as in the pp → ℓ + ℓ − γ analysis, if applicable.
In the √ s = 8 TeV analysis, both the photon transverse-momentum and the missing transverse-momentum cuts are increased to p γ T > 130 GeV and p νν T > 100 GeV, respectively. As in the Zγ → ℓ + ℓ − γ analysis the rapidity acceptance for jets is slightly increased to |η jet | < 4.5. The cuts are summarized in table 4.
The predicted cross sections at LO, NLO and NNLO can be found in table 5. The results are presented summed over all three neutrino species in the final state. In the inclusive case, i.e. for N jet ≥ 0, we find relatively large NLO corrections of around 57% and 68% and NNLO corrections of around 12% and 14% at √ s = 7 TeV and √ s = 8 TeV, respectively. The inclusive NNLO cross section prediction at √ s = 7 TeV is in good agreement with the cross section measured by ATLAS. In the exclusive case, N jet = 0, the NNLO corrections are very small, and the scale uncertainties are reduced down to the 1% level. We observe quite a significant discrepancy with respect to the ATLAS measurement for √ s = 7 TeV. Such discrepancy, however, is not completely unexpected. First of all, as mentioned in section 3.3 the stability of the fixed order calculation when a jet veto is applied is challenged and the perturbative uncertainties we find through scale variations are likely to be underestimated. Moreover, the Z → νν decay implies that the final state can be identified only through the photon and the additional radiation. Hadronization corrections, which are at the 1 − 2% level for all the other processes studied in this paper, lead in this case to sizeable effects for N jet = 0. The comparison of our NLO result with that quoted in table VII of ref.
[9], which is corrected for hadronization effects, indeed shows that in this case an O(30%) correction must be applied to the parton level theoretical prediction, thus reconciling it with the experimental result. Figure 7 shows the photon transverse-momentum and the missing transversemomentum distributions. These distributions are identical for Born kinematics due to momentum conservation, so the difference results purely from real-radiation corrections. Above the photon transverse-momentum cut of p γ T > 100 GeV, the difference between the two distributions is very small. Below a missing transverse momentum of 100 GeV, the cross section is only non-vanishing starting from the NLO. Figure 7 shows a perturbative instability around p T,miss ≈ 100 GeV. This instability originates from the incomplete cancellation of virtual and real corrections close to the phase space boundary (see ref. [84] for a discussion of this phenomenon.) This class of singularities is integrable and does not alter the inclusive cross section, but would require a resummed computation to achieve a reliable differential prediction close to the boundary.  Table 5. pp → ν ℓ ν ℓ γ cross sections at LO, NLO and NNLO. Scale uncertainties are evaluated as in table 2. The numerical uncertainty of the NNLO prediction from statistical error and finite r cut is conservatively estimated to be about 0.5%. The last column provides the corresponding result by ATLAS. We can also study the transverse-mass distribution of the ννγ system, defined as Figure 8 also shows the contribution from gluon fusion separately, which again is quite small and amounts to less than 2% of the fiducial cross section in the inclusive case and about 3% in the exclusive case.

pp → ℓν ℓ γ
We now present results for pp → ℓν ℓ γ at √ s = 7 TeV and 8 TeV. We again use the event selection criteria adopted in the ATLAS analysis [9]. This set of cuts is identical to that used in the pp → ℓ + ℓ − γ analysis, apart from the fact that the cut on the invariant mass of the leptons is replaced with a cut on the missing transverse momentum (which coincides with the transverse momentum of the neutrino from the W decay) of p ν T > 35 GeV. As in the case of Zγ in our √ s = 8 TeV analysis the rapidity of the jets is required to be |η jet | < 4.5. A summary of all cuts and event selection criteria can be found in table 6. All results in the following will be presented summed over the W charges, i.e. we combine the processes pp → W + γ and pp → W − γ, to facilitate the comparison with experimental data. The predicted fiducial cross sections both for the inclusive and the exclusive case can be found in table 7. In the inclusive case, the NLO corrections are quite large, and amount to about 136-143%. The NNLO corrections increase the NLO result by 19-20%. The impact of higher order corrections is thus much larger than in the case of Zγ production. We will come back to this point in section 3.6. Table 7 also shows the cross sections measured by ATLAS. The measurement of the inclusive cross sections shows a 2σ excess with respect to the NLO prediction, which is reduced to well below 1σ when including the NNLO corrections.
The impact of QCD corrections at NLO and NNLO is reduced to 60% and 7%, respectively, when the jet veto is applied (N jet = 0). As discussed in section 3.3, such an effect is expected and apparently leads to a more stable perturbative prediction, but also to the possible need of more conservative procedures to estimate perturbative uncertainties. In

JHEP07(2015)085
Frixione isolation with ε γ = 0.5, R = 0.4, n = 1 Jets anti-k T algorithm with D = 0.4  Table 7. W ± (→ ν ℓ ℓ) γ cross sections with the soft p γ T cut (p γ T > 15 GeV). Scale uncertainties are computed as in table 2. The numerical uncertainty of the NNLO prediction from statistical error and finite r cut is conservatively estimated to be about 0.8%. The last column provides the measured cross sections provided by ATLAS. the exclusive case, the excess of the measured fiducial cross sections over the theoretical prediction is reduced from 1.6σ to 1.2σ when going from NLO to NNLO. We note that the scale variations at NLO significantly underestimate the impact of the NNLO corrections, in particular in the inclusive case. Figure 9 shows the photon transverse-momentum distribution in comparison with the ATLAS measurement, both in the inclusive and in the exclusive case. Although the experimental uncertainties are large, the agreement between data and theory is clearly improved when including the NNLO corrections, in particular if no veto on jets is applied. Figure 10 shows the W γ cross section differential in the transverse mass of the ℓν ℓ γ system, normalized by the total fiducial cross section at the respective order.  Table 8. W ± (→ ν ℓ ℓ) γ cross sections with the hard p γ T cut (p γ T > 40 GeV). Scale uncertainties are computed as in table 2. The numerical uncertainty of the NNLO prediction from statistical error and finite r cut is conservatively estimated to be about 0.5%.
The calculation is done with the hard photon transverse-momentum cut p γ T > 40 GeV. The corresponding fiducial cross sections at LO, NLO, and NNLO are reported in table 8. The impact of QCD radiative corrections is 242-260% and 26% at NLO and NNLO, respectively. In figure 10, due to the normalization, the large overall corrections mostly cancel out, in particular at high transverse mass, and we observe only a slightly improved agreement with data when going from NLO to NNLO.
The increased relative impact of NLO and NNLO corrections when a harder p γ T cut (p γ T > 40 GeV) is applied can, in analogy to the Zγ case (see section 3.3), be better understood by studying the transverse-mass distributions with the soft and hard p γ T cut in more detail. The corresponding plots with a finer binning are shown in figure 11. When p γ T > 15 GeV, for Born kinematics the transverse mass has a lower bound at about m ℓνγ T 75 GeV, i.e. below the W → ℓν ℓ γ peak. When the photon transverse-momentum cut is increased to 40 GeV, the lower bound increases to m ℓνγ   Figure 11. Transverse-mass distribution of the ℓν ℓ γ system at LO (blue, dotted), NLO (red, dashed) and NNLO (green, solid) for p γ T > 15 GeV (left) andi p γ T > 40 GeV (right), in the inclusive case (N jet ≥ 0). The lower panel shows the NNLO/NLO ratio, and the bands indicate theoretical uncertainty estimates from scale variations. the region where the cross section is sizeable, and thus explaining the effect on the fiducial cross section.

The difference between W γ and Zγ
It is interesting to compare the relative size of the NLO and NNLO corrections to the Zγ and W γ processes we have considered. The results are summarized in table 9. It is clear that the W γ process features much larger radiative effects with respect to the Zγ processes. This should be contrasted to what happens in the case of inclusive W and Z boson production, where QCD radiative corrections are essentially identical [85]. It is thus the emission of the additional photon that breaks the similarity between the charged current and the neutral current processes. Table 9. Summary of the relative NLO and NNLO corrections in the channels under investigation,

JHEP07(2015)085
Numbers are reported for both soft and hard p γ T cuts.
Restricting the analysis to NLO for the moment, the main source for the difference between Zγ and W γ can be traced back to the gq and gq channels, which contribute a moderate, negative amount to the cross section in Zγ production, but are large and positive for W ± γ. It turns out that this effect is driven by resonant W γ contributions to the cross section, i.e. by pp → W (→ ℓν ℓ )γ topologies, and not by pp → W → ℓ(→ ℓγ)ν ℓ topologies, where the photon is emitted from the final-state lepton. These two contributions can only be separated in a gauge-invariant way if the W bosons are treated as on-shell particles, i.e. in a narrow-width approximation. By studying the LO contributions to the Zγ and W γ cross sections (see figure 1 and figure 2) it turns out that in W γ there is an additional Feynman diagram in which the photon is radiated off the W boson (see figure 2d). This additional diagram is responsible for a radiation zero [86], an exact zero present in the on-shell partonic W γ tree-level amplitude at cos θ * = 1/3, where θ * is the scattering angle in the centre-of-mass frame. This radiation zero gets diluted by the convolution with the parton densities and by off-shell effects, but it is responsible for the suppression of the Born level W γ cross section with respect to Zγ. Real radiation appearing at NLO breaks the radiation zero, and thus the relative impact of higher-order corrections is significantly increased. To quantitatively test this effect we consider the pp → ℓν ℓ γ and pp → ℓ + ℓ − γ processes studied in section 3.3 and section 3.5, with the same selection cuts. Contrary to what was done in the previous sections we disable the contributions from final state radiation and use the narrow-width approximation. In figure 12 (left) we plot the distribution in the rapidity difference ∆y ℓγ between the charged lepton and the photon [87].
We see that the LO distribution shows a pronounced dip at central rapidities. Although diluted by the convolution with the parton densities, the dip is clearly visible, and is responsible for the suppression of the W γ cross section. Since real radiation does not respect the radiation zero, the dip is filled up by radiative corrections. Roughly speaking, the NLO is a de facto LO prediction in the region of the dip and the NNLO corrections are thus relatively large as well. In contrast to W γ, the Zγ amplitude does not exhibit a radiation zero and, consequently no dip appears in the rapidity-difference distribution, as can be seen in figure 12 (right).
The presence of the radiation zero, and of the corresponding dip in the ∆y ℓγ distribution, are thus the reason for the increased importance of radiative corrections to the W γ processes.

Summary and discussion
In this paper we have presented the first complete and fully differential computation of QCD radiative corrections to W γ and Zγ production at hadron colliders. More precisely, we have considered the processes pp → ℓ + ℓ − γ, pp → ν ℓ ν ℓ γ and pp → ℓν ℓ γ, where, in the first case, the lepton pair ℓ + ℓ − is produced either by a Z boson or a virtual photon. The diagrams in which the photon is radiated off the final-state charged leptons were consistently included. We have presented quantitative predictions for fiducial cross sections for pp collisions at √ s = 7 and 8 TeV, and for various kinematical distributions (only at √ s = 7 TeV). The impact of QCD radiative corrections strongly depends on the applied cuts. In the case of

JHEP07(2015)085
Zγ, the impact of NNLO corrections is generally moderate, ranging from 8% to 18%. We have also shown that the loop induced gluon fusion contribution is generally small, and it accounts for less than 10% of the full O(α 2 S ) correction. In the case of W γ production the NNLO effects are more important, and range from 19% to 26%. The larger impact of QCD radiative effects in the case of W γ production is a well known consequence of a radiation zero [86] existing in the W γ amplitude at Born level. This effect produces a suppression of the LO distribution in the rapidity difference between the charged lepton and the photon, and NLO and NNLO corrections are thus quite significant. As expected, the impact of QCD radiative effects is strongly reduced when a jet veto is applied (N jet = 0), being smaller than 3% in the case of Zγ, and about 7% in the case of W γ.
We add few comments on the remaining uncertainties affecting our NNLO results. The uncertainties from missing higher-order contributions were estimated through scale variations, which are performed through independent variations of the renormalization and factorization scales around their central value (without constraints on their ratio). In the inclusive case the NNLO scale uncertainties obtained in this way are of the order of ±(1-2)% in the case of pp → ℓ + ℓ − γ (see table 2), ±(2-3)% in the case of pp → ν ℓ ν ℓ γ (see table 5), and ±4% in the case of pp → ℓν ℓ γ (see table 7). The comparison of the NNLO predictions to what is obtained at NLO shows that the NNLO-NLO difference is larger than the NLO scale dependence. We thus conclude that, as usual, scale variations can give only a lower limit on the true perturbative uncertainty. However, the NNLO is the first order at which all partonic channels are accounted for, and we believe that the NNLO scale uncertainties obtained in the case N jet ≥ 0 should provide the correct order of magnitude of the true uncertainty. As discussed, the situation is different for the case N jet = 0, in which the scale uncertainties are even smaller. Most likely, in this case a more conservative approach has to be adopted to obtain a realistic estimate of the perturbative uncertainty (see e.g. [80,83]).
The other source of uncertainty affecting our perturbative QCD calculations is the one coming from the PDFs. The PDF uncertainties at 68% CL that we obtain on our fiducial cross sections are at the 1%-2% level, both at NLO and NNLO, and are thus of the same order, or smaller, than the perturbative uncertainties. We have checked that, by using CT10 [88] and NNPDFs [89], the differences we obtain with the default MMHT result are of the same order.
The quantitative predictions we have presented for √ s = 7 TeV were obtained by using the same cuts adopted by the ATLAS collaboration in their measurement of the W γ and Zγ cross sections [9]. We have presented a comparison to ATLAS data, both for the fiducial cross sections and for some kinematical distributions. In the case of Zγ production the NNLO corrections slightly improve the agreement with the data, which, however, have large uncertainties. The only exception is the case pp → ν ℓ ν ℓ γ with a jet veto (N jet = 0), for which a 1.7σ discrepancy with the ATLAS result remains. In the case of W γ, the ATLAS result overshoots the NLO prediction by about 2σ. The large NNLO corrections reduce this excess to below 1σ. The NNLO corrections improve the agreement with the data also for the kinematical distributions we have studied, in particular for the p γ T distribution. However, it is known that this distribution, at large p γ T , is affected by JHEP07(2015)085 sizeable effects from EW corrections [23][24][25]. The impact of EW corrections depends on the way the additional photon radiation is treated. It is thus difficult at present to draw conclusions on the data/theory agreement in the high-p γ T region. More generally, we think it will be important to consistently combine the QCD calculations presented in this paper with a NLO computation of EW correction, which is, however, left for future work. [13] CMS collaboration, Measurement of the production cross section for Zγ → ννγ in pp collisions at √ s = 7 TeV and limits on ZZγ and Zγγ triple gauge boson couplings,