Advancing MiNNLO$_{\rm PS}$ to diboson processes: $Z\gamma$ production at NNLO+PS

We consider $Z \gamma$ production in hadronic collisions and present the first computation of next-to-next-to-leading order accurate predictions consistently matched to parton showers (NNLO+PS). Spin correlations, interferences and off-shell effects are included by calculating the full process $pp\to\ell^+\ell^-\gamma$. We extend the recently developed MiNNLO$_{\rm PS}$ method to genuine $2\to 2$ hard scattering processes at the LHC, which paves the way for NNLO+PS simulations of all diboson processes. This is the first $2\to2$ NNLO+PS calculation that does not require an a-posteriori multi-differential reweighting. We find that both NNLO corrections and matching to parton showers are crucial for an accurate simulation of the $Z \gamma$ process. Our predictions are in very good agreement with recent ATLAS data.


Introduction
Lacking clear hints for new-physics phenomena, particle phenomenology at the Large Hadron Collider (LHC) has entered the precision era. The accurate measurement of Standard Model (SM) processes provides a valuable alternative in the discovery of new-physics phenomena through small deviations from SM predictions. Many LHC reactions, in particular coloursinglet processes, are not only measured, but also predicted at a remarkable accuracy. For instance the recent Zγ [1] and ZZ [2] measurements include the full Run-2 data and are hitting percent-level uncertainties even for differential observables. The vast amount of data collected at the LHC will continuously decrease experimental uncertainties, thereby demanding accurate theory predictions in many relevant physics processes.
The theoretical description of fiducial cross sections and kinematic distributions has been greatly improved by the calculation of NNLO corrections in QCD perturbation theory. Those have become the standard for 2 → 1 and 2 → 2 colour-singlet processes  by now. The recent NNLO calculation of γγγ [36,37] production marks another milestone for precision calculations, since it is the first 2 → 3 genuine LHC process to be computed at this level of precision. A comparison of theory predictions to LHC data also highlights that a knowledge of NNLO corrections is crucial for theory results to describe data within their experimental uncertainties.
Vector-boson pair production processes in particular have become an integral part of the rich precision programme at the LHC. Being measured by reconstructing the vector bosons from their leptonic decay products, those processes offer clean experimental signatures with rather small experimental uncertainties. Apart from the measurement of their production rates and distributions, they provide a proxy for both direct and indirect searches for beyond-the-SM (BSM) physics. While very precise SM predictions of diboson processes are not needed to find a light resonance structure in invariant mass distributions, new-physics effects can give rise also to modifications and distortions of kinematic distributions. These effects, which can be parametrized through anomalous triple-gauge couplings or effective operators, enter vector-boson pair processes already at the leading order (LO). Constraining new-physics effects in these type of indirect searches crucially relies on accurate theory predictions for event rates and shapes of distributions.
Zγ production in the Z → + − decay channel provides a particularly pure experimental signature as the final state can be fully reconstructed. In combination with its relatively large cross sections this process is well suited for precision phenomenology. Indeed, Zγ production was measured extensively at the LHC at 7 TeV [38][39][40][41][42][43], 8 TeV [44][45][46][47], and 13 TeV [1,48], and ref. [1] was the first diboson analysis to include the full Run II data set. Even small deviations from the production rate or distributions in this process would be a direct hint of BSM physics. So far, full agreement with the SM was found, which provides a strong test of the gauge structure of electroweak (EW) interactions and the mechanism of EW symmetry breaking. On the other hand, the measurement of a non-zero ZZγ coupling, which is absent in the SM, would be direct evidence of physics beyond the SM (BSM). Moreover, Zγ final states are relevant in direct searches for BSM resonances and in Higgs boson measurements, see e.g. ref. [49,50], with the SM production being an irreducible background. Although the Higgs decay into a Zγ pair is rare, since it is loop-induced in the SM, effects from new-physics may significantly enhance this decay channel.
A substantial effort has been made in fixed-order calculations for Zγ production in the past years. Next-to-leading order (NLO) QCD corrections were computed some time ago both for on-shell Z bosons [51] and including their leptonic decays [52]. The first contribution known at next-to-next-to-leading order (NNLO) QCD was the loop induced gluon-fusion contribution [53][54][55]. Ref. [56] combined the NLO cross section, including photon radiation off the leptons, with the loop-induced gluon fusion contribution. The full NNLO QCD cross section for + − γ production was calculated in refs. [21,22] at the fully differential level and it was later confirmed in ref. [23] by an independent calculation. Also the NLO electroweak (EW) corrections are known [57,58].
However, the validity of fixed-order calculations is limited to observables dominated by hard QCD radiation. In kinematical regimes where soft and collinear QCD radiation becomes important, the perturbative expansion in the strong coupling constant is challenged by the appearance of large logarithmic contributions. The analytic resummation of those logarithms is usually restricted to a single observable, see e.g. the recent next-to-next-tonext-to-logarithmic (N 3 LL) results for the Zγ transverse momentum (p T, γ ) in refs. [59,60] for instance, or at most two observables, such as the joint resummation of logarithms in p T, γ and in the leading jet transverse momentum (p T,j 1 ) at next-to-next-to-logarithmic (NNLL) [61]. Parton showers, on the other hand, offer a numerical approach to include all-order effects in all phase-space observables simultaneously. Although their all-order logarithmic accuracy is rather limited (see ref. [62] for a recent discussion on this topic), parton-shower simulations are extremely important for experimental analysis since they allow for an exclusive description of the final state. Moreover, as measurements operate at the level of hadronic events, they require full-fledged parton-shower Monte Carlo simulations. In fact, any BSM analysis searching for small deviations from SM predictions at event level requires parton-shower predictions that include the highest possible fixed-order accuracy. While matching of NLO QCD predictions and parton showers (NLO+PS) has been worked out a while ago in seminal papers [63][64][65], 1 current experimental measurements demand the inclusion of NNLO QCD corrections in event generators to fully exploit LHC data.
So far four different NNLO+PS approaches have been presented in the literature [67][68][69][70][71], and all of them are formulated for colour singlet processes only. The methods of refs. [67,70,71] originate from the MiNLO procedure [67,72], which uses a transverse momentum resummation to upgrade an NLO calculation for X + 1 jet to become NLO accurate both for X and X + 1 jet. A numerical method which allows to extend the MiNLO procedure also to more complex final states was presented in ref. [73], where in particular the case of Higgs production in association with up to two jets was worked out. NNLO+PS approaches have been mostly applied to the simple 2 → 1 LHC processes, such as Higgs-boson production [70,71,74,75], the Drell-Yan process [69-71, 76, 77], Higgsstrahlung [78][79][80], which in terms of QCD corrections is still a 2 → 1 process, and to the 1 → 2 process H → bb decay [81].
Up to now, only the approach of ref. [67], which relies on a numerically highly demanding multi-dimensional reweighting in the Born phase space, has been applied to a genuine 2 → 2 process, namely W + W − production [82]. 2 This calculation has taken the reweighting procedure to its extreme. In fact, the Born phase space for W + W − → e + µ − ν eνµ involves nine variables (after taking the azimuthal symmetry into account). Ref. [82] had to recur to a number of features of the W -decays, such as the fact that the full angular dependence of boson decays can be parametrized through eight Collins-Soper functions, and used the fact that QCD corrections are largely independent of the off-shellness of the vector bosons to 1 For the processes considered in this work, in ref. [66] NLO QCD predictions for Z + γ + jets with different jet-multiplicities have been merged using the MEPS@NLO approach, which separates samples using a merging scale, and then interfaced to a parton shower. 2 The W + W − simulation is based on the MiNLO calculation of ref. [83], and the NNLO calculation of ref. [30] performed within the Matrix framework [84].
reduce the number of kinematic variables in the parametrization of the Born phase space. Still, the residual Born-phase space dependence had to be discretized. The finite bin sizes used in the reweighting limits the accuracy of the results in regions of phase space characterized by coarse binning. This is typically the case in regions close to the kinematic edges and limits of the phase space, such as tails of kinematic distributions, which, on the other hand, are particularly interesting for BSM searches. Actually, the numerical limitations related to an a-posteriori reweighting constitutes a problem already for a way simpler process such as Drell-Yan, since in this case, given the very large amount of data available, experiments require a high theoretical precision over the whole phase space.
In this paper, we consider Zγ production and present the first NNLO+PS calculation for a genuine 2 → 2 process that includes NNLO corrections directly during event generation, without any post-processing or reweighting of the events being required. In fact, this is also the first time a NNLO Zγ calculation independent of a slicing cutoff is performed (cf. refs. [21][22][23]). To this end, we have extended the just recently developed MiNNLO PS method [70] to deal with genuine 2 → 2 reactions, which paves the way for NNLO+PS simulations of all other diboson processes. As anticipated already in ref. [70], MiNNLO PS is a very powerful approach as its underlying idea applies beyond 2 → 1 processes and beyond colour-singlet production, with the following features: • NNLO corrections are calculated directly during the generation of the events, with no need for further reweighting.
• No merging scale is required to separate different multiplicities in the generated event samples.
• When combined with transverse-momentum ordered parton showers, the matching preserves the leading logarithmic structure of the shower simulation. 3 We consider all topologies leading to the final state + − γ in our calculation, including offshell effects and spin correlations. Since, neither a Zγ nor a Zγ+jet generator was available within the POWHEG-BOX framework [85], we also present new calculations of these two processes using the POWHEG method [64]. Our implementation of the Zγ+jet generator builds the basis for the inclusion of NNLO QCD corrections to Zγ production through the MiNNLO PS method. The ensuing calculation allows us to retain NNLO QCD accuracy in the event generation and to interface it to a parton shower, which is a necessary step for a complete and realistic event simulation. In particular, multiple photon emissions through the QED shower, as well non-perturbative QCD effects using hadronization and underlying event models can be included. While these corrections have only very mild effects on fully inclusive quantities, they can have a substantial impact on jet-binned cross-sections and other more exclusive observables measured at the LHC. This manuscript is organized as follows: in section 2 we discuss the implementation of Zγ and Zγ+jet generators within the POWHEG-BOX-RES framework, where we introduce the two processes (section 2.1), some basics about POWHEG-BOX-RES (section 2.2), and give details about the treatment of photon isolation and its practical implementation (section 2.3). The inclusion of NNLO corrections to Zγ production into Zγ+jet generator is discussed in section 3, including details on the extraction of the the two-loop amplitude (section 3.1), a general discussion on how we extended the MiNNLO PS approach to 2 → 2 processes (section 3.2), and further practical details relevant for the specific case of the Zγ MiNNLO PS generator (section 3.3). In section 4, after describing the setup used in our calculation and the set of fiducial cuts used in the analysis (section 4.1), we first present some validation plots (section 4.3), compare to the most accurate predictions for the transverse momentum distribution of the colourless system (section 4.4), and finally compare our results for fiducial cross sections and distributions to ATLAS data (section 4.5). We conclude and summarize in section 5. Some technical aspects are discussed in more details in the appendices.

NLO+PS simulation of Zγ and Zγ+jet production
In this section, we discuss the implementation of NLO+PS generators for Zγ and Zγ+jet production in the POWHEG-BOX framework [85]. Both processes were not yet available in this framework and we present their first calculation in the POWHEG [65] approach at NLO+PS. Moreover, the Zγ+jet process serves as starting point to reach NNLO accuracy for Zγ production through the MiNNLO PS method, as detailed in section 3. Since these processes involve an EW resonance, we exploit the POWHEG-BOX-RES code [86], which is specifically designed to deal with intermediate resonances. In the following, we first introduce the two processes under consideration, then we recall some relevant features of the POWHEG-BOX-RES framework, and finally discuss details regarding the treatment of the photon in the final state and the QED singularities associated to it.

Description of the processes
We consider the production processes where ∈ {e, µ} is a massless charged lepton. For brevity, we refer to these processes as Zγ and Zγ+jet production in the following. As illustrated in figure 1, Zγ production is initiated by quark-antiquark annihilation at LO. The photon can be emitted either by the quark line (q-type diagrams) or by the lepton line ( -type diagrams), each of which yields a different resonance structure of the respective amplitudes. Sample LO diagrams for Zγ+jet production are shown in figure 2, with the same classification into q-type and -type diagrams. The distinction between those two resonance structures will be relevant when treating them as two different resonance histories Figure 1: Sample LO diagrams for the + − γ production with two different resonance stuctures. Figure 2: Sample LO diagrams for + − γ+jet production including q-type diagrams (a-c) and -type diagrams (d-e).
within the POWHEG-BOX-RES framework, discussed in section 2.2. In addition to the tree-level Born amplitudes, the NLO calculation of the Zγ (Zγ+jet) process requires the respective one-loop contributions as well as the tree-level real emission Zγ+jet (Zγ+2-jet) amplitudes.
The NLO corrections to Zγ and Zγ+jet production have been implemented within the POWHEG-BOX-RES framework. For the Zγ generator the relevant amplitudes have been extracted from MCFM [87], while for the Zγ+jet generator those have been implemented both using MCFM and via an interface to OpenLoops 2 [88]. The helicity amplitudes of MCFM are implemented from the analytic expressions computed in refs. [89,90] for Zγ production and in ref. [91] for Zγ+jet production. The width of the Z boson is included in the fixed-width scheme. For Zγ+jet, the contribution from third generation quarks inside loops has been entirely removed for those diagrams where the Z boson attaches to a fermion loop through an axial-vector coupling, while the massless bottom loop has been retained for those contributions where the corresponding top effects decouple as 1/m top 4 [23]. The impact of this approximation is expected to be rather small as shown in ref. [92], where the leading heavy-quark loop contribution has been evaluated in the 1/m top 2 expansion in the context of Zj and Zjj. We further note that, in view of the NNLO calculation for Zγ production discussed in section 3, omitting the contribution of third generation quarks is in line with the fact that the heavy-quark loop contributions at two loops are currently not known, and therefore not included throughout our NNLO+PS results.
When using the OpenLoops interface for the Zγ+jet computation of the NLO amplitudes, the complex mass scheme can be used and the full top-mass effects can be accounted for, while in the MCFM amplitudes the width is implemented only in a fixed-width scheme and heavy-quark loop effects are included only approximately. Since QED effects are included just at LO, the difference between the complex-mass and the fixed-width scheme amounts to an overall normalization, whose impact is below 0.1%. When comparing results with full top-mass effects as available in Openloops to approximate results as implemented in MCFM we found per mille effects for quantities inclusive on QCD radiation. This is expected, since heavy-quark effects at one loop are non-vanishing only in the presence of final-state radiation. For jet-related quantities, the discrepancy between the two approaches can range from a few percent at low transverse momentum up to approximately 10-20% in highly boosted regions (p T > 2m top ) or high-invariant mass regions. This is not surprising, since the process at hand involves s-channel fermion-loop contributions, which become more important in these phase space regions. For observables involving a jet our results are NLO accurate only, hence characterized by larger theoretical uncertainties. In summary, we find that mass effects are always much smaller than our quoted theoretical uncertainties and for the numerical studies performed in this paper, which are not devoted to boosted regions, it is perfectly fine to use approximate results for the heavy-quark mass effects. Accordingly, because of the better numerical performance of MCFM, we choose to use these amplitudes to obtain the results of this paper. Specifically, we find that MCFM virtual amplitudes are about ten times faster than Openloop ones. On the other hand, one can make use of the folding option in POWHEG, where the real contribution is evaluated multiple times for each virtual one. This improves the numerical performance whenever the virtual amplitudes constitute the bottleneck in the numerical evaluation. For greater flexibility, in the release of the numerical code, the option to choose between the OpenLoops and the MCFM implementations has been made available.

The POWHEG-BOX-RES framework
We calculate NLO+PS predictions for Zγ and Zγ+jet using the POWHEG method, which is based on the following master formula [64,65,85]: where Φ B is the Born phase space of the process under consideration. The functionB(Φ B ) describes the inclusive NLO process, where extra QCD emissions are integrated out; the content of the curly brackets is responsible for the exclusive generation of one extra QCD radiation with respect to the Born process according to the POWHEG method [64,65,85]. B and R denote the squared tree-level matrix elements for the process at Born level and for its real radiation, respectively. The evaluation of the POWHEG Sudakov form factor ∆ pwg [64] depends on the real phase space Φ rad through the transverse momentum p T,rad of the extra radiation. The POWHEG cutoff Λ pwg is used to veto QCD emissions in the non-perturbative regime and its default value is Λ pwg = 0.89 GeV. The parton shower then adds additional radiation to eq. (2.2) that contributes beyond NLO with respect to the underlying Born at all orders in perturbation theory. We refer to the original publications for explicit fomulae [64,65,85].
For the practical implementation of the Zγ and Zγ+jet generators we exploit the POWHEG-BOX-RES framework [86], which takes into account the different resonance structures of each process. All possible resonance histories are automatically identified and the code performs a resonance-aware phase space sampling. The efficiency of the infrared subtraction is improved by means of its resonance-aware subtraction algorithm, where the mapping from a real to the underlying Born configuration preserves the virtuality of intermediate resonances [86,93]. Moreover, in a parton-shower context the distortion of resonances through recoil effects is avoided by supplying it with details on the resonance cascade chain. The resonance awareness of the phase space sampling and the subtraction improve the numerical stability of the evaluation of theB function in the POWHEG formula of eq. (2.2).
The key idea behind the algorithm used in the POWHEG-BOX-RES framework is to decompose the cross section into contributions associated to a well-defined resonance structure, which are enhanced on that particular cascade chain. As discussed in section 2.1, both Zγ and Zγ+jet production have two different resonance histories, which can be associated to q-type diagrams, where the photon is emitted from the quark/antiquark line, and -type ones, where the final state photon is radiated off one of the two leptons.
The POWHEG-BOX-RES framework takes as input bare flavour structures B of the Born process (e.g. B = {uū → + − γ}), which only contain the information on the initial- After introducing the full flavour and resonance structureˆ B of the Born process (e.g.
, which embodies details on the entire resonance history, we can further decompose B B as a weighted sum overˆ B using weight functions Pˆ where T ( B ) is named a tree and denotes all graphs with the given initial-and final-state flavour configuration B . The weight functions Pˆ B are chosen such that eq. (2.4) expresses B B as a sum over resonance-peaked terms Pˆ B Bˆ B , which develop the expected resonance enhancement ofˆ B . There is a certain freedom in their explicit expression, and in the POWHEG-BOX-RES code the following choice is made: where the sum in the denominator runs over all configurations in the tree T ( B (ˆ B )) of graphs having a bare flavour structure B (ˆ B ) compatible withˆ B . In the specific case of Zγ and Zγ+jet production, Pˆ B assumes two different functional forms depending on whether B refers to q-type or -type diagrams, and they read where s is the invariant mass of the lepton pair and s γ that of the produced colour-singlet system. The same discussion applies to the virtual corrections, while for the real-emission contribution a similar decomposition is performed taking into account the different singular regions. As described in detail in ref. [86], the concept of resonance history directly affects the definition of QCD singular regions: only soft/collinear singular regions compatible with a given full real flavour structureˆ R are considered in the FKS decomposition of the real amplitude. It is important to note that, for each of the full real flavour structureŝ R the mappings from the real to the Born configurations preserve the virtualities of the intermediate resonances, which is crucial to guarantee a cancellation of singularities between real corrections and their counterterms.

Treatment of the isolated photon and details of the implementation
The emission of a soft or collinear photon from a quark or a charged lepton induces QED singularities. Processes with final-state photons therefore require not only suitable criteria to isolate photons in the experimental analyses, but they also call for an IR-safe isolation procedure on the theory side. Since in the POWHEG-BOX framework fiducial cuts are usually applied at analyses level after parton showering, which modifies the kinematics of the final states, we discuss how to include photon-isolation requirements already at the event generation level in this framework to obtain IR-safe predictions.
To produce isolated photons in the final state there are two relevant mechanisms: the direct production in the hard process, which can be described perturbatively, and the production through fragmentation of a quark or a gluon, which is non-perturbative. The separation between the two production mechanisms in theoretical predictions is quite delicate, as sharply isolating the photon from the partons would spoil infrared (IR) safety. So-called fragmentation functions are required to absorb singularities related to collinear photon emissions in the latter production mechanism. Those functions are determined from data with relatively large uncertainties. On the other hand, Frixione's smooth-cone isolation of the photons [94] offers an IR-safe alternative that completely removes the fragmentation component. This substantially simplifies theoretical calculations of processes with isolated photons at higher orders in perturbation theory. Although the direct usage of smoothcone isolation in experimental analyses is not possible due to the finite granularity of the calorimeter, data-theory comparisons are facilitated by tuning the smooth-cone parameters to mimic the fixed-cone isolation used by the experiments, see e.g. ref. [95].
So far, only few processes involving final-state photons have been implemented in the POWHEG-BOX framework: W γ [96] and the direct photon [93,97]. These two generators make use of the photon fragmentation component. In particular in ref. [96] the hadron fragmentation into photons is modelled within POWHEG in combination with a QCD+QED shower. In this case the theoretical prediction can apply directly the photon isolation criteria imposed by the experiments to distinguish prompt photons taking part in the hard scattering process from possible background sources (such as photons from decay of π 0 mesons or from quark fragmentation). This facilitates a direct comparison between experimental and theoretical results. Those isolation criteria limit the hadronic activity in the vicinity of the photon by imposing where the sum of the transverse energy E had T of the hadrons inside a fixed cone of radius R 0 around the photon is constrained to be less than E max T .
In view of the NNLO extension considered in this paper, presented in section 3, we instead rely on smooth-cone isolation [94] to turn off the fragmentation component and to deal with QED collinear singularities perturbatively in an IR-safe manner. In this case, phase-space configurations where the photon becomes collinear to a quark are removed while preserving IR safety by allowing arbitrarily soft QCD radiation within smoothly decreasing cones around the photon direction. In practice, this means that the smooth-cone isolation is implemented by restricting the hadronic (partonic) activity within every cone of radius δ = (∆η) 2 + (∆φ) 2 < δ 0 around a final-state photon, where δ 0 sets the maximal cone size, by imposing the following condition such that the total hadronic (partonic) transverse energy inside the cone is smaller than E max T (δ). The parameter n controls the smoothness of the isolation function and E ref T is a reference transverse-momentum scale that can be chosen to be either a fraction γ of the transverse momentum of the respective photon (p T,γ ) or a fixed value (p 0 T ), In our calculations we impose smooth-cone isolation on the phase space of all Zγ+jet and Zγ+2-jet configurations. Furthermore, various technical phase-space cuts at event generation level are necessary in order to obtain IR safe results. Those generation cuts and parameters of the smooth-cone isolation are given in appendix A. They are chosen to be much looser than the ones eventually applied at analyses level after parton showering. We stress that, since we also employ suppression factors for the NLO squared amplitudes (as discussed in detail below), the resulting differential cross section times suppression factors vanishes in the singular regions, which will not pass fiducial cuts.
As commonly used in many POWHEG-BOX generators, we exploit the possibility to split the real squared matrix element R into a singular and a finite (remnant) contribution. Such splitting improves the numerical performance of the code, especially as far as the efficiency of the event generation is concerned, in cases where the ratio R/B departs from its corresponding soft/collinear approximation, for instance in presence of Born zeros [98]. Following section 5 of ref. [85], we write the splitting into a singular and a remnant contribution for each singular region α of the real amplitude as: where Φ R is the real phase-space and S is called damping factor. Only the singular contribution α R α sing. (Φ R ) is exponentiated in the POWHEG Sudakov ∆ pwg and used to generate the first emission according to the POWHEG method in eq. (2.2), while the finite remnant contribution α R α remn. (Φ R ) can be treated separately, by generating it with standard techniques and feeding it directly into the parton shower.
A standard damping factor [85] is used in both Zγ and Zγ+jet generators, where S=0 when the real squared amplitude in a singular region is greater than five times its soft/collinear approximation, and S=1 otherwise. Additionally, to improve the numerical convergence, the Zγ+jet generator requires a special setting of the damping factor, which ensures that QED singularities appearing in the real squared matrix element are moved into the remnant contribution. Indeed, not all of the QED singular regions appearing in the real matrix elements have a singularity in their underlying Born. Accordingly, the associated real singularity is not mitigated by an overall Born suppression factor (as described in more detail below). To deal with this issue, we define, in each singular region, the invariant mass m rad of the emitter-emitted pair of that singular region, which is the quantity that becomes small close to QCD singularities, and we use as a damping factor 4 where c ∈ [0, 1] is a free parameter that we choose below, and the sum runs over all (initialand final-state) quarks with when i is a quark in the initial state , The splitting induced by the suppression factor in eq. (2.11) is such that, when a QED singularity dominates, the event is included in R α remn. (Φ R ) (S → 0), and, when the QCD singularity is dominant, the event is moved into R α sing. (Φ R ) (S → 1). The numerical constant c controls the transition region between QED singularities in d iγ and QCD singularities in m 2 rad . Since α e /α s ∼ 0.1, we use the value c = 0.1, which we have checked to be suitable for an efficient generation of events.
Finally, we exploit another tool of the POWHEG-BOX framework that allows us to improve the numerical convergence in the relevant phase-space regions. By introducing suppression factors, which multiply the cross section during integration and are a posteriori divided out again, it is possible to redirect the numerical sampling of events into certain regions in phase space. This is mandatory for processes that have singularities at Born level, such as the QED singularities in Zγ and Zγ+jet production, to avoid sampling large statistics in phase-space regions which are eventually removed by the fiducial cuts at analysis level. In order to obtain suitable integration grids that give more weight to the phase-space regions selected by the fiducial cuts, we have introduced a Born suppression factor that vanishes in singular regions related both to QCD and QED emissions. Its precise form is given in appendix A. Since the real phase space is generated directly from the Born one in the POWHEG-BOX, the same factor is also applied to R α sing. (Φ R ). For the remnant contribution R α remn. (Φ R ), on the other hand, which is QCD regular, but is in our case QED singular, an analogous suppression factor has been chosen, that is given in appendix A as well.
We stress that the implementational details discussed throughout this section are by no means standard, despite the fact that the tools we are using existed already in the POWHEG-BOX framework. The proper adjustment of those tools and their related parameters required a great effort, which was necessary to obtain a generator for Zγ+jet production with sufficient numerical efficiency to be extended to NNLO corrections to Zγ production discussed in the next section.
3 Reaching NNLO accuracy for Zγ production using MiNNLO PS We use our implementation of the Zγ+jet generator in the POWHEG-BOX-RES framework discussed in section 2 as a starting point. To include NNLO corrections for Zγ production in this calculation we employ the recently developed MiNNLO PS method [70]. To this end, we extend the MiNNLO PS method to processes with non-trivial two-loop corrections, i.e. genuine 2 → 2 hard-scattering processes such as vector-boson pair production. After a general discussing of the ingredients required for a NNLO calculation, we recall the details of the MiNNLO PS method and present its extension to 2 → 2 processes. Finally, we provide practical details on how that calculation is embedded in the POWHEG-BOX-RES framework.

Ingredients of a NNLO calculation
In section 2.1 we have discussed the contributions relevant to evaluate NLO corrections to Zγ and Zγ+jet production. Those involve tree-level and one-loop amplitudes for the processes pp → Zγ and pp → Zγ+jet as well as the tree-level amplitude for pp → Zγ+2jet. The same amplitudes enter the NNLO calculation for Zγ production, i.e. at the Born level and as real, virtual one-loop, real-virtual and double-real corrections. The only missing ingredients for the NNLO calculation are the two-loop corrections in the qq channel (sample diagrams are shown in figure 3 (a-b)), and the loop-induced contribution in the gluon-fusion channel (with a sample diagram shown in figure 3 (c)). The latter can be separated from the others since it is effectively only LO accurate. Its size is rather small, being less than 10% of the NNLO corrections and below 1% of the full Zγ cross section at NNLO [84]. We thus refrain from including the loop-induced gluon-fusion contribution in our calculation. We note, however, that while its calculation at LO+PS is quite straightforward and easily done with current standard tools, see e.g. refs. [85,99], a more sophisticated treatment would require to match the NLO corrections to the loop-induced gluon-fusion contribution with parton showers. Despite being feasible with current technology, this is beyond the scope of this paper and left for future studies.
For the two-loop corrections we use the qq → γ helicity amplitudes calculated in ref. [100]. Those have been fully implemented into the Matrix framework [84,101] using the results of ref. [100]. In order to exploit this implementation and evaluate the two-loop helicity amplitudes within our MiNNLO PS calculation presented in this paper we have compiled Matrix as a C++ library and linked it to our POWHEG-BOX-RES Zγ+jet generator.
We further exploit Matrix for all fixed-order NNLO results used for comparisons throughout this paper. The Matrix framework features NNLO QCD corrections to a large number of colour-singlet processes at hadron colliders. Several state-of-the-art NNLO QCD predictions [21,22,25,26,[28][29][30][31][32]37] 5 have been obtained with this framework, and for massive diboson processes it has been extended to combine NNLO QCD both with NLO EW corrections [103] and with NLO QCD corrections to the loop-induced gluon fusion contribution [104,105]. Through the recently implemented Matrix+RadISH interface [59,61] the code now also combines NNLO calculations with high-accuracy resummation through the RadISH formalism [106][107][108] for different transverse observables, such as the transverse momentum of the colour singlet.
3.2 Generalization of MiNNLO PS to 2 → 2 colour-singlet processes In the following we recall the central aspects of the MiNNLO PS method of ref. [70] as well as some of its improvements presented in ref. [71], and we discuss in detail all the changes required to apply the method to a genuine 2 → 2 colour-singlet process, such as vectorboson pair production. The main difference compared to 2 → 1 hadronic processes, such as Higgs and Drell-Yan production, is that the one-and two-loop virtual corrections for a general process can not be written as a simple form factor times the Born amplitude, and thereby receive a dependence on the respective flavour configuration. To this end, we will recast, where appropriate, the relevant formulae of refs. [70,71] in a flavour-dependent notation.
MiNNLO PS is a method to perform a NNLO calculation for a produced colour singlet F with invariant mass Q fully differential in the respective Born phase space Φ F , in such a way that it can be subsequently matched to a parton shower. In the context of this paper, F would be the Zγ (or more precisely + − γ) final state, but we prefer to keep the discussion general throughout this section. We employ the same notation as in refs. [70,71] in the following. The starting point of MiNNLO PS is a POWHEG implementation of colour singlet production in association with one jet (FJ), whose phase space we denote by Φ FJ . We thus write eq. (2.2) explicitly with the Born process being the FJ production: Here,B(Φ FJ ) describes the FJ process, i.e. including the first radiation, using the full NLO cross section while the content of the curly brackets accounts for the second QCD emission through the POWHEG mechanism, with B and R being the squared tree-level matrix elements for FJ and FJJ production, respectively, and Φ rad (p T,rad ) referring to the phase space (transverse momentum) of the second emission. Radiation beyond the second one is generated by the parton shower, which adds corrections of O(α 3 s (Q)) and higher at all orders in perturbation theory. In order to reach NNLO accuracy in the phase space of the color singlet F in eq. (3.1), we modify the content of theB(Φ FJ ) function, which is the central ingredient of the MiNNLO PS method.
The derivation of theB(Φ FJ ) function in MiNNLO PS [70] is based on the following formula that describes the transverse momentum of the color singlet (p T ) up to NNLO and refs. [70,71,80,82].
is fully differential in the Born phase space Φ F : where R f includes only non-singular contributions at small p T , and At variance with the formulas in refs. [70,71], we have introduced an explicit sum over the flavour structures F of the Born process pp → F. The quantities without index F should be understood as having been summed implicitly over F as in the original formulation of refs. [70,71], in particular The flavour dependence of these quantities originates entirely from the hard-virtual coefficient function H F , 7 which contains the virtual corrections that, for a general 2 → 2 hadronic process, become dependent on the flavour and on the Born phase-space Φ F . Up to two loops it is given by The convolution between two functions f (z) and g(z) is defined as Note thatH F in eq. (3.5) just includes an additional shift with respect to H F , see eq. (4.26) of ref. [70].
This dependence propagates through theB (2) F coefficient to the Sudakov form factor, since the derivation of the MiNNLO PS formalism is based on setting the renormalization scale µ R ∼ p T , which exponentiates the H (1) F coefficient and requires a redefinition ofB (2) F , as discussed in the derivation of the replacement in eq. (4.26) of ref. [70]: . A few comments are in order: all quantities with index F also depend on the Born kinematics. For ease of notation we do not indicate explicitly their Φ F dependence. The G functions [109] in eq. (3.5) are present only in the case of gluon-induced reactions, i.e. they are zero for Zγ production and kept here only for completeness. For a colour-singlet process the Born flavours F correspond to the two initial-state partons, which have been denoted with cc in eq. (3.5). In that equation |M F | 2 cc denotes the Born matrix element squared,C are the collinear coefficient functions, and f F is given in section 3.3. Writing also the differential NLO cross section for FJ production as a sum over its where [X] (i) denotes the coefficient of the i th term in the perturbative expansion of the quantity X, allows us to recast our starting formula in eq. (3.2) as where F ← FJ denotes the projection from the flavour structure of FJ production to the one of the Born process F. This projection is trivial in the case of Zγ production, as the Born is always qq initiated and only the respective quark-flavour is of relevance in the F dependence ofS F (p T ). Furthermore, we have introduced a new symbol D F (p T ) that accounts for the relevant NNLO corrections with the following two choices of treating terms beyond accuracy. In the original MiNNLO PS formulation of ref. [70] we have truncated eq. (3.11) to third order in α s (p T ), i.e.
With this truncation at the differential level eq. (3.11) does not reproduce anymore the exact total derivative that we started with in eq. (3.2). In order to preserve the total derivative and keep into account additional terms beyond accuracy, a new prescription has been suggested in ref. [71] which will be our default choice throughout this paper. The relevant expressions for its evaluation, including the ones of the [D F (p T )] (i) coefficients, are reported in appendix C and D of ref. [70] and in appendix A of ref. [71], where the flavour dependence can be simply included through the replacements F . Eq. (3.11), upon integration over p T , is NNLO accurate with both choices of D F (p T ), as they differ only by terms of O(α 4 s ) and higher. As discussed in detail in refs. [70,71], this is achieved by the consistent inclusion of all singular terms to third order in α s (p T ).
We can now make the connection to POWHEG and the parton-shower matching formula in eq. (3.1). The formulation of eq. (3.11) applies also to the fully differential cross section in the Φ FJ phase space, and it can be used to achieve NNLO accuracy for theB(Φ FJ ) function [70,71] where F corr FJ (Φ FJ ) guarantees a proper spreading of the Born-like NNLO corrections in the full Φ FJ phase space, as discussed in detail in Section 3 of ref. [70] and briefly recalled in the next section.

Practical details of the implementation within POWHEG-BOX-RES
We have applied the MiNNLO PS formalism discussed in the previous section to our implementation of the Zγ+jet generator in POWHEG-BOX-RES. In the following we discuss practical aspects of our calculation in that framework and set F = Zγ from now on.
We briefly recall how the NNLO corrections, which have Born kinematics and additionally depend on p T , are included in Zγ+jet generator. The relevant kinematics is obtained by performing a phase-space projection Φ ZγJ → Φ Zγ and by determining p T from the given Φ ZγJ phase-space point. Furthermore, the Born-like NNLO corrections need to be distributed in the Φ ZγJ kinematics. There is a certain degree of freedom in how to associate (Φ Zγ , p T ) with the full Φ ZγJ phase space. This is encoded in the factor F corr (Φ ZγJ ) of eq. (3.14) , (3.15) and in the details of the ZγJ → Zγ projection. The functions J ZγJ (Φ FJ ) for ZγJ = {qq, qg, gq} have been chosen according to the collinear limit of the tree-level matrix element squared of the Zγ+jet process, using eq. (A.14) of ref. [70]. This is a suitable compromise between computational speed and physically sound results, as the spreading is performed according to the pseudorapidity distribution of the radiation described by that approximation. With this choice the numerical convergence of the integral in the denominator of eq. (3.15) does not scale with the complexity of the process, which is a crucial requirement for multi-leg processes, such as the one at hand. We stress that the spreading factor in eq. (3.15) is designed in such a way that the integral over the Φ ZγJ phase space of its product with any function of the Φ Zγ Born variables yields exactly the integral of that function, when integrated over the Φ Zγ phase space (see eq. (3.2) of ref. [70]). Details of our projection to Zγ events are given in appendix B. Note that this projection does not preserve the full Born kinematics. While it keeps all invariant masses and the rapidity of the Zγ system unchanged, it does alter for instance the transverse momentum of the photon. As a result, the photon transverse momentum after the Φ ZγJ → Φ Zγ projection is neither bounded from below by the technical generation cuts nor controlled by the phase-space suppression factor introduced for the Φ ZγJ kinematics in section 2.3. This induces a singular behaviour through the Born and virtual amplitudes in both the Sudakov form factor and the luminosity factor defined in eqs. (3.5) and (3.6). We have therefore added the requirement p T,γ ≥ 10 GeV in the projected Φ Zγ kinematics. This technical cut is below the p T,γ threshold used at analysis level and it can be controlled through the input card. Its effect is strictly beyond accuracy, affecting only regular contributions at large p T, γ . In fact, as discussed in appendix B, our projection can lead to configurations with p T,γ → 0 GeV only for events where the jet is back to back to the Zγ system, and the Z and the photon are aligned with each other. It is then clear that for such events p T,j > p T,γ and there are no large logarithms associated to p T,j . Indeed, we have varied the cutoff down by a factor of ten, finding changes at the level of the numerical precision of less than a 1%. Embedding the MiNNLO PS corrections within the POWHEG-BOX-RES framework requires some further discussion in respect to the resonance-aware features. Since thē B(Φ Zγ ) function is modified in an additive way in eq. (3.14), such that the NNLO corrections are treated on the same footing as all other contributions with Φ ZγJ kinematics, we decompose them as a weighted sum overˆ B using the weight functions Pˆ B just like in eq. (2.4). Our final formula for the Zγ MiNNLO PS generator reads

(3.16)
In the following we provide some details on the treatment of higher-order terms and the scale settings. The discussion summarizes briefly the one in section 3.2 of ref. [71], which we refer to for more details. In the large transverse-momentum region it is important to switch off all-order resummation effects to avoid spurious contributions. As is standard, we do it by introducing modified logarithms via the replacement where we set here p = 6. As pointed out in ref. [71], in order to preserve the total derivative of eq. (3.2), this prescription requires three further adjustments. The lower integration bound of the Sudakov is to be replaced by p T → Qe −L ; D F (p T ) (or D F (p T ) (3) ) needs to be multiplied by a proper jacobian factor, see eq. (13) and (16) of ref. [71]; the perturbative scales need to be set consistently with the scale of the modified logarithms at small transverse momenta. Additionally, we smoothly approach non-perturbative scales at small p T and introduce the non-perturbative parameter Q 0 to regularize the Landau singularity, setting the central renormalization and factorization scales to [71] where we set Q 0 = 0.5 GeV, and g(p T ) has been chosen so as to suppress the Q 0 shift at large values of p T : The scale setting of eq. (3.18) provides a consistent treatment in the small p T region and preserves the total derivative of eq. (3.2). However, at large p T it yields µ R,0 ∼ µ F,0 ∼ Q, while a scale setting of µ R,0 ∼ µ F,0 ∼ p T might be preferred in that region. To this end, the scales entering the NLO Zγ+jet cross section in eq. (3.16) can be set via the flag largeptscales 1 to which we apply in our calculation. It is important that eq. (3.20) matches the scales of the Sudakov form factor and the D F (p T ) terms at small p T . At the same time it ensures a dynamical scale choice of µ R,0 = µ F,0 ∼ p T at large p T . Next, we specify the precise definition of the first and second order hard-virtual function of eq. (3.8) for Zγ production in our subtraction scheme. First, let us recall that all tree-level and one-loop amplitudes have been extracted from MCFM and for comparison also linked through OpenLoops, and that the (one-loop and) two-loop qq → γ amplitudes have been obtained by creating an interface to their implementation in Matrix. At variance with the 2 → 1 processes considered in ref. [70], it is not possible to provide compact expressions for the hard function of Zγ production. For brevity, we thus start from the expressions of the first and second order hard-virtual function where IR singularities have been subtracted according to the q T -scheme (more precisely, choosing the hard-scheme [110] as resummation scheme) that we denote as H qT(1)  in the following. Those coefficients are unambiguously defined in eqs. (12) and (62) of ref. [110] and can be extracted from the one-loop and two-loop virtual amplitudes using the expressions of that paper. In fact, the Matrix inteface directly provides the hard-virtual coefficients in that scheme, so that we only need to perform the appropriate scheme conversion to match the MiNNLO PS conventions [111]: 8 where C (1),δ cc and C (2),δ cc are the terms proportional to δ(1 − z) of the first and second order coefficients of the collinear coefficient function, which in the case of a quark-induced process (cc = qq) are given bỹ where C F = 4/3 and C A = 3, and N f is the number of light quark flavours. Note that C (1/2) qq in the MiNNLO PS convention can be obtained from the ones of ref. [110] by simply addingC (1/2),δ qq × δ(1 − z). For completeness, we also provide the corresponding expressions for the gluon-induced case (cc = gg) Finally, we conclude this section by reporting two further non-standard settings related to the showering of the Zγ MiNNLO PS events. First, we turn on by default the POWHEG-BOX doublefsr option, which was introduced and discussed in detail in ref. [112]. When this option is turned on, the emitter-emitted role is exchanged in events for which the emitter would be softer than the emitted particle. This considerably reduces the appearance of spikes in distributions due to events with large weights that pass fiducial cuts after showering. Furthermore, for the Pythia8 shower [113], we set the flag SpaceShower:dipoleRecoil 1 (see ref. [114]). Its effect is to have a global recoil (i.e. a recoil which affects all final state particles including the color-singlet system) only for initialinitial colour-dipole emissions, while for initial-final ones a local recoil scheme is used. The reason for this choice is that, as shown in ref. [71], a local recoil for initial-final emissions reduces the effect of the shower on the colour-singlet kinematics, in particular in large rapidity regions.

Phenomenological results
In this section, we present NNLO+PS accurate predictions for Zγ production. We consider different fiducial selections discussing both integrated cross sections and differential distributions in presence of fiducial cuts. MiNNLO PS predictions are compared to NNLO and MiNLO results. This allows us to validate the accuracy of our predictions for observables inclusive over QCD radiation and observables requiring the presence of jets. At the same time, we demonstrate the importance of both NNLO accuracy in the event simulation and the inclusion of parton-shower effects. Finally, we compare MiNNLO PS predictions to a recent ATLAS measurement [1].

Input parameters, settings and fiducial cuts
We present predictions for 13 TeV collisions at the LHC. Our results have been obtained by using the G µ -scheme, where the electroweak coupling is defined as where cos 2 θ W = m W 2 /m Z 2 and the input parameters are set to We use N f = 5 massless quark flavours, and we choose the corresponding NNLO PDF set of NNPDF3.0 [115] with a strong coupling constant of α s (m Z ) = 0.118. In the case of the fixed-order predictions we use the PDF set at the respective order in QCD perturbation theory. To be precise, for MiNLO and MiNNLO PS we read the PDF grids using the lhapdf interface [116], copy them into hoppet grids [117] and then use the hoppet evolution code to perfom the DGLAP evolution of the PDFs, also for scales below the internal PDF infrared cutoff, as pointed out in section 3.3. The calculation of D F (p T ) in eq. (3.13) requires the evaluation of different PDF convolutions and the computation of polylogarithms. For the latter we made use of the hplog package [118]. Our setting of the central renormalization and factorization scales (µ R,0 , µ F,0 ), which is in line with the MiNNLO PS (MiNLO ) method, has been discussed at length in section 3.3, see eqs. (3.18) and (3.20). In all fixed-order results presented throughout this section we adopt the following setting of the central renormalization and factorization scales: where m γ is the invariant mass of the Zγ system. The different scale choice between a NLO/NNLO fixed-order and a MiNLO /MiNNLO PS calculation induces effects beyond the nominal accuracy, in addition to the different treatment of higher-order terms, see comments close to eq. (3.13). As a results, minor differences between the fixed-order and matched predictions are expected even for more inclusive cross-sections. Nevertheless, the results should largely agree within scale uncertainties, at least in cases where scale uncertainties are expected to be a reliable estimate of missing higher-order corrections.
We employ 7-point scale variations as an estimate of the theoretical uncertainties. Accordingly, the minimum and maximum of the cross section have been taken over a set of scales (µ R , µ F ) obtained by varying K R and K F in µ R = K R µ R,0 and µ F = K F µ F,0 , respectively, within the values: (2,1), (1,2), (1, 1), (1, 1/2) , (1/2, 1) , (1/2, 1/2)} , (4.4) We reiterate that, when performing scale variations for MiNLO and MiNNLO PS , we also include the scale dependence of the Sudakov form factor to account for additional sources of uncertainties [70], which are absent in a fixed-order calculation. For all technical details and choices made for our implementation of the Zγ+jet generator as well as for the treatment of the NNLO corrections through the MiNNLO PS method we refer the reader to section 2.3 and in section 3.3, respectively. In particular, the settings discussed in section 2.3 are essential to get a good convergence of the Monte Carlo integration and an efficient event generation, by adopting generation cuts and individual suppression factors at Born level, for the singular real contributions and the remnant contributions (cf. also appendix A). In addition, the folding of the radiation variables (ξ, y, φ) [85,119] has been used (with a choice of (1, 5, 1) for the folding parameters) to evaluate the double-real correction (Zγ+2-jet) more often, which further improves the numerical convergence. Despite all those efforts to achieve a better numerical performance, we had to produce ∼ 100 million Zγ events with our MiNNLO PS generator to obtain acceptable statistical uncertainties and predict integrated cross sections in the fiducial setups considered here at the level of a few permille. Still, our comparison of differential distributions to NNLO predictions suffers from some fluctuations. We note, however, that with Zγ production we have picked probably the most involved diboson processes, featuring various complications, in particular its substantial complexity with respect to the QED singularity structure. We therefore expect other diboson processes to have a much lower demand for numerical resources.
As pointed out before, we omit the loop-induced gluon-gluon component in our implementation and throughout this paper. This contribution is relatively small, being only ∼ 1% of the NNLO cross section, and it can be incoherently added to our predictions through a dedicated calculation, which, however, is beyond the scope of this paper. Consequently, we drop the loop-induced gluon fusion contribution also from the fixed-order calculation, when comparing our MiNNLO PS against NNLO results, to warrant a meaningful comparison.
We employ the Pythia8 parton shower [113] with the one of the A14 tunes [120] (specifically py8tune 21) to dress the hard event with further soft/collinear QCD radiation with the default POWHEG setting for the parton-shower starting scale. The showered results do not include any effects from hadronization or underlying event models. Moreover, the photon is required to be generated only at the hard scattering level: contributions from a QED shower or the decay of unstable particles is not included. Finally, we keep photons stable by preventing any photon conversion effect, i.e. no γ → + − or γ →qq splitting.
We present results for two sets of fiducial cuts, which are summarized in table 1. The first one is identical to that used in refs. [22,84] and motivated by in an earlier ATLAS analysis [43]. We refer to it as ATLAS-setup-1 in the following. The second one was instead ATLAS-setup-1 [43] ATLAS-setup-2 [1] Lepton cuts p T, > 25 GeV |η | < 2.47 Frixione isolation with Frixione isolation with n = 1 γ = 0.5 δ 0 = 0.4 used in the most recent ATLAS 13 TeV measurement of ref. [1] using the full Run-2 data and named ATLAS-setup-2 in the following. We make use of ATLAS-setup-1 for validation purposes and to show the importance of NNLO+PS matching, while ATLAS-setup-2 is also used to compare MiNNLO PS predictions with the most updated experimental measurement available for Zγ production. Both setups in table 1 involve standard transverse momentum and pseudorapidity thresholds to identify leptons and photons as well as a lower invariant-mass cut on the lepton pair. ATLAS-setup-2 places an additional requirement on the sum of the invariant masses of the Zγ system and of the lepton pair. This cut suppresses the contribution from -type diagrams, where the photon is emitted from the final state leptons (cf. figure 1b), enhancing t-channel production through q-type diagrams (cf. figure 1a). Moreover, separation cuts between two particles (i, j) are applied in ∆R ij = ∆η 2 ij + ∆φ 2 ij , where ∆η ij and ∆φ ij are their differences in the pseudorapidity and the azimuthal angle, respectively. In both setups leptons are separated from the isolated photon, while only ATLAS-setup-1 imposes an additional separation of jets from leptons and from the isolated photon, which in turn requires a jet definition. As a consequence, we employ ATLAS-setup-1 to study jet observables and show NLO/LO accuracy of MiNNLO PS predictions for Zγ+jet/Zγ+2-jet configurations. Finally, isolation criteria for the photon are needed, as detailed in section 2.3, which is done by means of Frixione isolation in both setups. In ATLAS-setup-2, Frixione isolation in a smaller cone and a second isolation criterium is applied by requiring the scalar sum of the transverse energy of all stable particles (except neutrinos and muons) within a cone around the photon of size R = 0.2 to be less than 7% of the photon transverse momentum (see [1] for more details). Note that we apply the latter isolation criterium only when analyzing events after parton showering, but not at Les-Houches-Event (LHE) or fixed-order level.

Fiducial cross sections
In table 2 we report predictions for the Zγ cross section in the two fiducial setups at LO, NLO and NNLO, and for MiNLO and MiNNLO PS matched to Pythia8. The fixed-order results have been obtained with Matrix. Although MiNNLO PS and NNLO calculations entail a different treatment of terms beyond accuracy, in both setups the agreement of their predicted cross sections is remarkably good. One should bear in mind, however, that in ATLAS-setup-2 there is a slight difference in the treatment of the isolated photon at fixed order, which does not include the E cone0.2 T /p T,γ < 0.07 cut, as discussed in the previous section. We further notice from table 2 that the scale uncertainties of the theoretical predictions are increased in ATLAS-setup-2. This is caused by the additional m + m γ cut and the larger cut on the photon transverse momentum in that fiducial setup, rendering the predictions more sensitive to additional QCD radiation that is described at a lower perturbative accuracy.
Comparing MiNNLO PS and MiNLO predictions, the inclusion of NNLO corrections through MiNNLO PS has a rather moderate effect for the fiducial cross section of +2% in ATLAS-setup-1 and +0.5% in ATLAS-setup-2. In fact, in both cases (and particularly evident for the latter setup) MiNLO predictions are actually closer to the NNLO results than to the NLO ones. After all, the Sudakov form factor in eq. (3.6) is exactly the same for MiNNLO PS and MiNLO , and MiNLO predictions already contain various contributions beyond NLO accuracy, including all real corrections at NNLO through the merging of NLO corrections to Zγ+jet production. Still, by reaching NNLO accuracy through MiNNLO PS the predictions get even closer to the NNLO results and the uncertainty bands substantially decrease, by almost a factor ten in ATLAS-setup-1 and by more than a factor 2 in ATLAS-setup-2. Indeed, the MiNNLO PS scale uncertainties are comparable with the NNLO ones. The fact that they are slightly larger is expected since the MiNNLO PS pro-cedure probes lower scales both in the PDFs and in α s , and it includes scale variations also for the Sudakov form factor.
Finally, we find a very good agreement of our NNLO+PS accurate MiNNLO PS predictions with the cross section measured by ATLAS in ref. [1], which are perfectly compatible within the quoted experimental errors. We now turn to discussing differential distributions in the fiducial phase space. In this section we compare our MiNNLO PS predictions with MiNLO and NNLO results. This serves two purposes. On the one hand, MiNNLO PS distributions are validated for one-jet and two-jet observables against the ones obtained with MiNLO and for Born-level observables (inclusive over QCD radiation) against NNLO predictions. On the other hand, this allows us to show the importance of NNLO+PS matching with respect to less accurate results. To these ends, we discuss selected distributions which are particularly significant to show the performance of MiNNLO PS predictions. The figures of this and the upcoming sections are organized as follows: the main frame shows the predictions from MiNNLO PS (blue, solid line), together with all other results relevant for the given comparison. In an inset we display the bin-by-bin ratio of all the histograms that appear in the main frame to the MiNNLO PS one. The bands indicate the theoretical uncertainties that are computed from scale variations.
We start by discussing quantities that involve jets in the final state in figure 4 and figure 5, where MiNNLO PS (blue, solid line) and MiNLO (black, dotted line) are compared at the un-showered LHE level in ATLAS-setup-1. Since for these observables MiNLO and MiNNLO PS have the same accuracy the two predictions are expected not to differ from each other significantly (i.e. not beyond uncertainties), both in terms of shapes and size of scale uncertainty bands. In particular, such agreement serves as a validation that NNLO corrections are properly spread by the factor in eq. (3.15) in the jet-resolved phase space of Zγ+jet production without altering the NLO accuracy. As a matter of fact, the left plot of figure 4 shows that MiNNLO PS and MiNLO predictions agree well within uncertainties for both the pseudorapidity difference between the Zγ system and the hardest jet (∆η γ, j 1 ). Furthermore, the size of the uncertainty bands are comparable over the whole pseudorapidity range. In a similar manner, the ratio between MiNNLO PS and MiNLO is nearly flat for the invariant mass of the photon and the hardest jet (m γj 1 ) in the right plot of figure 4. Here, we further observe the effect of the Frixione isolation, which dampens the distribution in the photon-jet collinear limit. Also in figure 5 MiNNLO PS and MiNLO predictions agree well for the distance between the photon and the leading and subleading jet in the η-φ plane (∆R γj 1 and ∆R γj 2 ). As ∆R γj 2 involves the second-hardest jet, both MiNNLO PS and MiNLO are only LO accurate, which is also evident from the broadening of the uncertainty bands. We have examined a large number of other quantities involving jets (not shown here) observing a similar behaviour in all cases. Next, in figure 6, figure 7 and figure 8, we compare MiNNLO PS against MiNLO (black, dotted) and NNLO predictions from Matrix (red, dashed) for Born-level observables (inclusive over QCD radiation) after showering with Pythia8. By and large, we observe a very good agreement of MiNNLO PS and NNLO predictions, especially considering the fact that they differ from each other in the choice of the renormalization and factorization scales, and in the treatment of higher-order contributions. What can be appreciated is the clear reduction of the scale uncertainties of MiNNLO PS predictions with respect to the MiNLO ones up to a size which is comparable to the NNLO ones.
In particular, figure 6 displays the pseudorapidity distribution of the Zγ system (η γ ) in each of the two fiducial setups. The ratio of NNLO over MiNNLO PS is close to one in both cases, with uncertainty bands of one to two percent in ATLAS-setup-1. In ATLAS-setup-2, on the other hand, the bands are roughly twice as large, as already observed for the integrated cross section due to the higher sensitivity to phase-space regions related to real QCD radiation.
In figure 7 we show the distributions in the invariant mass (m ) and transverse momentum (p T, ) of the lepton pair in ATLAS-setup-1. We can appreciate the Z-boson resonance in the m distribution as well as a broader, but smaller, enhancement around m ∼ 70 GeV, caused by the Z-boson resonance in m γ . The qualitative behaviour of MiNNLO PS with respect to MiNLO and NNLO predictions in the ratio inset are relatively similar to the one of the Zγ rapidity distribution. However, there seems to be a small effect induced by the resummation of the parton shower reducing the enhancement due to m γ , where the NNLO prediction is slightly larger. In the p T, distribution on the right we observe an interesting behaviour of the NNLO prediction. The NNLO result develops a perturbative instability (Sudakov shoulder) [121] around p T, ∼ 15 GeV caused by an incomplete cancellation of virtual and real contributions from soft gluons, which is logarithmically divergent, but integrable. The reason is the fiducial cut p T,γ > 15 GeV (see table 1) that for LO kinematics implies p T, = p T,γ > 15 GeV, so that the p T, distribution is not filled below 15 GeV at LO.  Thus, the fixed-order result is NNLO accurate only for p T, > p T,γ , while for p T, < p T,γ at least one QCD emission is necessary, which is described only at NLO accuracy. At the same time, at threshold the prediction becomes sensitive to soft-gluon effects, resulting in an instability at fixed order. Indeed, the parton shower cures this behaviour and yields a physical prediction at threshold for both MiNLO and MiNNLO PS . This is one example where an NNLO calculation is insufficient and NNLO+PS matching is required.
In figure 8 we consider distributions in the second-hardest lepton, showing its transverse momentum (p T, 2 ) in the left and its rapidity (η 2 ) in the right plot. Similar conclusions as made before for m and η γ apply also for these observables, so no further comments are needed. We reiterate however that, while the central predictions of MiNLO and MiNNLO PS are generally close to each other, since MiNLO already includes many terms beyond NLO accuracy for Zγ production, scale uncertainties are substantially reduced in case of MiNNLO PS down to the level of the NNLO ones. We continue our discussion of differential distributions with the transverse-momentum spectrum of the Zγ system (p T, γ ). In figure 9 we compare the p T, γ distribution in ATLAS-setup-2 obtained with MiNNLO PS against a more accurate prediction at NNLO+N 3 LL (green, dashed curve), using the analytic resummation of large logarithmic contributions within Matrix+RadISH [59,61]. For comparison we also show the NNLO result, which is effectively only NLO accurate for this distribution. The NNLO+N 3 LL prediction uses [59] µ R,0 = µ F,0 = m 2 + p T,γ 2 and Q res,0 = 1 2 m γ (4.5) as central scales, where Q res,0 is the central resummation scale, which is varied by a factor of two up and down, while taking the envelope together with the 7-point µ R and µ F variation for the total scale uncertainty. The p T, γ spectrum is shown in two different ranges in figure 9. From the wider range in the left plot we notice that despite the different scale settings in the three calculations their predictions are in reasonable agreement at large p T, γ values. In this region all predictions are only effectively NLO accurate, which is indicated by the enlarged scaleuncertainty bands. At small p T, γ the fixed-order result becomes unphysical, as the distribution is logarithmically divergent in the p T → 0 limit, which is visible already in the left plot of figure 9, but can be appreciated morein the zoomed version on the right. In this region, only the calculations that properly account for the resummation of soft QCD radiation by means of an analytic procedure (NNLO+N 3 LL) or through parton shower simulations (MiNNLO PS ) provide a meaningful description. Even though at small p T, γ the Matrix+RadISH computation is N 3 LL accurate, while the parton shower has a lower logarithmic accuracy, MiNNLO PS and NNLO+N 3 LL predictions are in excellent agreement down to the transverse-momentum values (almost) in the non-perturbative regime.

Comparison of differential distributions against ATLAS data
Finally, we employ our MiNNLO PS generator to compare NNLO+PS accurate predictions directly to ATLAS results from the recent 13 TeV measurement of ref. [1], which relies on the full 139 fb −1 Run-2 data. The comparison, carried out in ATLAS-setup-2, is presented in figure 10. The experimental data are given as green points, with error bars that refer to the experimental uncertainty. Six observables are shown: the transverse momentum (p T,γ ) and the pseudorapidity (η γ ) of the photon, the transverse momentum (p T, γ ) and the invariant mass (m γ ) of the Zγ system, together with their ratio p T, γ /m γ and the difference in the azimuthal angle between the lepton pair and the photon (∆φ , γ ).
Overall, we observe a remarkably good agreement both in the predicted shapes of the distributions and in the normalization, especially given the fact that the theoretical and the experimental uncertainties are at the few-percent level only. All data points agree with our predictions within the experimental error bars, with the exception of only very few bins, where the agreement is reached within twice the experimental error. This is a clear improvement over the NLO-accurate event simulations employed in the data-theory comparison in figure 6 of ref. [1], both in terms of accuracy (i.e. to describe the data) and in terms of precision (i.e. regarding theoretical uncertainties). Moreover, looking at the comparison of NNLO-accurate predictions to data in figure 7 of ref. [1], it is clear that some (more inclusive) observables are equally well described at fixed order, while for observables sensitive to QCD radiation, such as p T, γ and ∆φ , γ , NNLO predictions are not sufficient, and the matching to a parton shower is essential. In conclusion, our MiNNLO PS calculation combines the two most important aspects (NNLO and parton-shower effects) to provide the most accurate and most precise Zγ predictions to date, which will be essential to find potential deviations from the SM for this process in future.
Let us discuss in more detail the p T, γ distribution in figure 10. In this plot we have also added the more-accurate NNLO+N 3    the scale setting of eq. (4.5) that induces differences also at large p T, γ . Despite the good agreement of MiNNLO PS with data, the analytically resummed result is performing even better, especially in the first few bins, where the higher accuracy in the resummation of large logarithmic contributions is important. Although MiNNLO PS and NNLO+N 3 LL agree quite well (cf. the discussion in section 4.4), this shows that for an observable like p T, γ it can be very useful to resort to tools that predict a single distribution more accurately, if available. Nevertheless, it is reassuring that our accurate multi-purpose MiNNLO PS simulation, with all its flexibility to predict essentially any IR-safe observable, provides a very good description of such distributions as well.
We further notice that the deviation at small m γ is due to missing QED effects as shown in ref. [1]. Our MiNNLO PS computation renders the inclusion of such effects with a NNLO prediction feasible by using a QED shower within Pythia8, which could be very useful in an experimental analysis. This is however beyond the scope of this paper and left for future studies.

Summary
We have presented a novel calculation of NNLO+PS accurate predictions to Zγ production at the LHC. This is the first calculation of a genuine 2 → 2 process at this accuracy that does not require an a-posteriori multi-differential reweighting. To this end, we have extended the MiNNLO PS method [70] to 2 → 2 colour-singlet production, with non-trivial one-loop and two-loop corrections, and applied it to the Zγ process. More precisely, we have considered all resonant and non-resonant contributions to the hard-scattering process pp → + − γ, including off-shell effects and spin correlations.
As a starting point we have implemented NLO+PS generators for both Zγ and Zγ+jet production within the POWHEG-BOX-RES framework [86]. The Zγ+jet generator has then been extended to include NNLO corrections to Zγ production by means of the MiNNLO PS method. The two-loop virtual corrections [100] are linked through their implementation within the Matrix code [84], which we treat as a library. A substantial amount of work has been devoted to render those calculations numerically efficient, and we have discussed the technical and physical aspects in detail. All our computations will be made publicly available within POWHEG-BOX-RES. 9 We have presented NNLO+PS accurate predictions using our MiNNLO PS calculation in hadronic collisions at 13 TeV. Observables exclusive in the final state jets have been used to validate the way we include NNLO corrections through the MiNNLO PS procedure by comparing MiNNLO PS with MiNLO predictions. Observables inclusive over QCD radiation are generally well described by a fixed-order NNLO calculation. Although MiNNLO PS predictions differ in the treatment of terms beyond accuracy, we found them to be in very good agreement with NNLO results, both for fiducial cross sections and for differential distributions. Moreover, we have shown the importance of NNLO+PS accurate predictions. On the one hand, they render predictions physical for observables sensitive to soft-gluon effects, where NNLO results fail to provide a suitable description. On the other hand, the inclusion of NNLO corrections through MiNNLO PS achieves a substantial improvement over MiNLO results in terms of scale uncertainties for inclusive observables. We also compared MiNNLO PS predictions for the Zγ transverse-momentum spectrum against a more accurate analytically resummed calculation at NNLO+N 3 LL and found a remarkable agreement down to transverse-momentum values close to the non-perturbative regime. Finally, we have shown that MiNNLO PS predictions are in excellent agreement with the latest ATLAS 13 TeV data, with various improvements over lower-order simulations.
This calculation paves the way for NNLO+PS predictions for all diboson processes in the future. Moreover, the contribution from the loop-induced gluon fusion component, despite being rather small for Zγ production, could be calculated separately at (N)LO+PS and added to our predictions. A rather straightforward advancement would also be to consider Z → νν decays in future, which will be highly relevant also for dark-matter searches. In this context, also suitable modifications of the SM could be considered, for instance by introducing effective Z Zγ and γ γZ couplings.
Finally, we reckon that the presented results as well as our MiNNLO PS generator for Zγ production will be a useful advancement over previous Monte Carlo predictions and tools. Especially since the MiNNLO PS generator can be directly applied in experimental analyses of Zγ production or searches, the improved theoretical predictions in terms of accuracy and precision might stimulate further studies of this process in future.

A Generation cuts and suppression factors
Since our Born process involves a number of QED and QCD singularities, we make use of Born and remnant suppression factors to sample the phase space. Additionally, we introduce a number of small technical cuts in the phase space generation. In this Appendix we give all details about the generation cuts and suppression factors that we have used to obtain the results presented in this paper.
We start by outlining the generation cuts that we employ. First, we introduce a lower cut p cut T,γ = 5 GeV on the photon transverse momentum, which is required to avoid QED singularities related to collinear photon emissions from the initial state. We also impose a similar cut of p cut T,j = 1 GeV on the transverse momentum of the outgoing QCD parton. A lower cut m cut = 40 GeV on the invariant mass of the lepton pair is imposed to avoid singular configurations in γ * → splittings. Note that, since the invariant mass of the resonances are preserved when radiation is generated, any value m cut equal or below the cut used in the analysis is allowed. Furthermore, we require that the photon is isolated from leptons and QCD partons in the final state. For this purpose, we introduce a cut m 2 lγ = 0.1 GeV 2 , where l = {e + , e − } and we introduce a smooth isolation as in eq. (2.8) with E ref T = γ p T,γ with δ 0 = 0.05, γ = 0.5 and n = 1. All these generation cuts can be modified via the input card, provided their values is much smaller than the values used in the calculation of the fiducial cross sections.
The Born suppression factor that we adopt is constructed in factorized form with F supp (p T,γ ) = (p T,γ ) 2 (p T,γ ) 2 + (p cut T,γ ) 2 , with p cut Since we apply a MiNLO procedure, we do not need any suppression related to the Born outgoing parton. It is clear that, whenever a singularity is approched, the Born suppression factor in eq. (A.1) will vanish in such a way that the cross section times Born suppression factor itself will remain finite. As discussed in section 2.3, we have a remnant contribution which is QCD regular but QED singular. Accordingly, we introduce a remnant suppression factor of the form R supp = F supp (p T,γ ) · G supp (∆R γ,l + ) · G supp (∆R γ,l − ) · H supp (∆R γ,j ) (A.5) · H supp (∆R γ,j2 ) · H supp (∆R j1,j2 ) · L supp (p t,j2 )· , (A.6) with L supp (p t,j2 ) = (p t,j2 ) 2 (p t,j2 ) 2 + (p cut T,j ) 2 , with p cut T,j = 20 GeV . (A.7) As usual in POWHEG, the Born suppression is evaluated using the Born kinematics, while the remnant suppression is evaluated using the real radiation partonic kinematics. We note that in our case it is not necessary to introduce the additional factor L supp (p t,j2 ), however we found that results converged more quickly with the introduction of this additional suppression factor.

B Projection to the Zγ underlying Born configuration
The evaluation of the last term in Eq. (3.16) requires a projection from the Zγ+jet to Zγ kinematics. Here we give details about this projection and comment on configurations, which, after projection, have p T,γ close to zero.
We denote by p 1 and p 2 the two incoming momenta, and by p γ , p Z and p j the momenta of the photon, Z boson and jet in the final state. We define p tot = p 1 + p 2 − p j = p γ + p Z . Our projection to the underlying Born configuration consists of a longitudinal boost(by β L ), such that, after boosting, p tot has no z component. Then, a second boost (by β T ) in the transverse plane, such that p tot after both boosts has no transverse component, is applied, followed by a final boost back in the longitudinal direction (by −β L ). We add a prime to all quantities after the first longitudinal boost and a double prime to those after the second one. The boost vector of the transverse boost is then given by After the second boost, the transverse momentum of the photon becomes p T,γ = p T,γ γ T p T,γ + β T E γ , with γ T = 1 Therefore, after this boost, the condition p T,γ = 0 is met if By comparing Eq. (B.1) and (B.3) we see that p T,γ and p T,j must be anti-aligned. Furthermore, since E tot = E γ + E Z > E γ we have that p T,j > p T,γ . Accordingly, any boost that leads to a vanishing transverse momentum of the photon in the underlying Born configuration has a jet that is harder than p T,γ in the ZγJ configuration. Since for the photon we impose a transverse momemtum cut, this region is free of any large logarithm, and corrrections such as D F (p T ) in eqs. (3.13) and (3.14) are beyond our accuracy and can be dropped.