MiNNLO$_{\text{PS}}$: A new method to match NNLO QCD to parton showers

We present a novel method to combine QCD calculations at next-to-next-to-leading order (NNLO) with parton shower (PS) simulations, that can be applied to the production of heavy systems in hadronic collisions, such as colour singlets or a $t\bar{t}$ pair. The NNLO corrections are included by connecting the MiNLO$^\prime$ method with transverse-momentum resummation, and they are calculated at generation time without any additional reweighting, making the algorithm considerably efficient. Moreover, the combination of different jet multiplicities does not require any unphysical merging scale, and the matching preserves the structure of the leading logarithmic corrections of the Monte Carlo simulation for parton showers ordered in transverse momentum. We present proof-of-concept applications to hadronic Higgs production and the Drell-Yan process at the LHC.


Introduction
Particle phenomenology at the Large Hadron Collider (LHC) has entered the precision era. After the landmark discovery of the Higgs boson [1,2], which explains the electroweak (EW) symmetry breaking and completes the particle content predicted by the Standard Model (SM), the LHC is now focussing upon the search for hints of new-physics phenomena. Despite several indications that there must be physics beyond the SM (BSM) at relatively low scales, as of now new-physics searches at the LHC have been unsuccessful.
In this scenario, precision measurements have become of foremost importance, as they enhance the sensitivity of indirect searches for new physics through small deviations from SM predictions.
Modifications of the electroweak sector, induced by many BSM theories, can be unravelled through precision studies of electroweak interactions. Both the production of a Higgs-boson and of EW vector bosons play a crucial role in this respect. Notably, for the Drell-Yan process the experimental uncertainties have already surpassed the percent level [3][4][5]. On the theory side, predictions have to be controlled at the same level of precision, which demands calculations with at least next-to-next-to-leading order (NNLO) accuracy in QCD perturbation theory. Furthermore, the experimental measurements operate at the level of hadronic events and require full-fledged Monte Carlo simulations in their analyses. The inclusion of NNLO QCD corrections in event generators is therefore mandatory to fully exploit LHC data.
In this paper, we present a novel method to perform the consistent matching of NNLO calculations and parton showers (hereafter NNLO+PS), based on the structure of transversemomentum resummation. Our method builds upon MiNLO [6], a procedure for improving NLO multijet calculations with the appropriate choice of scales and with the inclusion of Sudakov form factors, that is particularly suited to be interfaced with parton-shower generators using the POWHEG method. In ref. [7] the MiNLO procedure was refined in such a way that, in processes involving the production of a massive colour singlet system in association with one jet, the NLO accuracy is formally retained also for observables inclusive in the jet. Such procedure, dubbed MiNLO , yields an NLO multi-jet merging method that does not use a merging scale. Notably, a numerical method to extend the MiNLO procedure to more complex processes, and its application to Higgs production in association with up to two jets, was presented in ref. [8].
In this article we extend the MiNLO method to achieve NNLO accuracy at the fully differential level in the zero-jet phase space, while retaining NLO accuracy for the one-jet configurations. Our method will be referred to as MiNNLO PS in the following, and it has 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.
Maintaining the logarithmic accuracy of the shower is a crucial requirement of all NLO+PS, and a fortiori NNLO+PS approaches. We stress that this requirement is immediately met by the MiNNLO PS approach, that works by generating the first two hardest emissions and letting the shower generate all the remaining ones. We also recall that, if the shower ordering variable differs from the NLO+PS one, maintaining the Leading Logarithmic accuracy of the shower becomes a delicate issue. An example is given by a POWHEG based generator interfaced to an angular-ordered shower. To preserve the accuracy of the shower, not only one needs to veto shower radiation that has relative transverse momentum greater than the one generated by POWHEG, but also one has to resort to truncated showers to compensate for missing collinear-soft radiation. Failing to do so spoils the shower accuracy at leadinglogarithmic level (in fact, at the double-logarithmic level). 1 Three different NNLO+PS approaches have been previously formulated in the literature [7,14,15] and applied to the simplest LHC processes, namely Higgs-boson production [16,17] and the Drell-Yan process [15,18,19]. The approach of ref. [7] shares all features listed above except the first one, i.e. it requires a multi-dimensional reweighting of the MiNLO' samples in the Born phase space to achieve NNLO accuracy. It has been recently applied to more complicated LHC processes, such as the two Higgs-strahlung reactions [20,21], and the production of two opposite-charge leptons and two neutrinos (W + W − ) [22]. 2 These computations have employed the reweighting procedure to its extreme. Despite yielding physically sound results, the reweighting in the high-dimensional Born phase space of these processes poses substantial technical limitations. Apart from the numerical demand of the reweighting itself and certain approximations that had to be made in these calculations, the discretisation of the Born phase space through finite bin sizes of the reweighted observables reduces the applicability of the results in phase-space regions with coarse binning, usually located in the least populated regions of phase space (e.g. in the tails of the kinematic distributions). In fact, the numerical limitation of the reweighting constitutes a problem already for the simpler Drell-Yan process, since the experiments require a considerably large number of generated events for the current and future LHC analyses.
The MiNNLO PS method presented in this paper lifts these shortcomings, while retaining the same advantages of a MiNLO computation. The terms relevant to achieve NNLO accuracy are obtained by connecting the MiNLO formula with the momentumspace resummation of the transverse-momentum spectrum formulated in refs. [26,27]. This allows us to make a direct link between the resummation and the POWHEG procedure [9], resulting in a consistent NNLO+PS formulation.
Due to the substantially improved numerical efficiency compared to the reweighting approach, the MiNNLO PS method allows us to tackle without any approximations the full class of complex color-singlet final states, such as the highly relevant four-lepton (vectorboson pair production) processes. As proof-of-concept applications, we consider hadronic Higgs production and the Drell-Yan process at the LHC, and compare our results against previous predictions. These computations are implemented and will be made publicly available within the POWHEG-BOX framework [9,28,29]. 3 All-order, higher-twist, and nonperturbative QCD effects are modelled through the interface to a parton shower generator 1 Truncated shower were first discussed in ref. [9], and first implemented in the Herwig++ context in ref. [10]. Currently Herwig++ implements them in its internal NLO+PS, POWHEG-scheme processes (see ref. [11]). They have also been used in a slightly different context in refs. [12,13]. 2 The W + W − simulation is based on the MiNLO calculation of ref. [23], and the NNLO calculation of ref. [24] performed within the Matrix framework [25]. 3 Instructions to download the code will be soon made available at http://powhegbox.mib.infn.it.
which provides a realistic simulation of hadronic events. Despite the fact that the formulae presented here are limited to the hadro-production of heavy colour-singlet systems, our formalism is quite general, and can be applied to other processes, such as the production of heavy quarks.
The manuscript is organized as follows: In section 2 we describe in general terms the main idea behind the MiNNLO PS approach, and determine the relevant corrections to the MiNLO formulation necessary to reach NNLO accuracy. Practical aspects of the implementation of these new terms within the MiNLO framework are discussed in section 3. In section 4 we provide a more rigorous derivation of the MiNNLO PS method by starting from the momentum-space resummation formula for the transverse-momentum spectrum. Our proof-of-concept computations for Higgs and Z-boson production are presented in section 5, where we provide a full validation against existing results. We summarize our findings in section 6. A number of technical details and explicit formulae are summarized in appendices A to E.

Description of the procedure
In this section we describe the procedure to perform a consistent matching of a NNLO QCD calculation for the production of a heavy colour-singlet system to a fully exclusive parton-shower simulation. We start by recalling the necessary elements of the MiNLO method in section 2.1 and 2.2, while in section 2.3 we derive the additional terms necessary to achieve NNLO accuracy.

The MiNLO method
We review now the basic elements of the MiNLO method, and how it achieves NLO accuracy. We formulate it in a way that is as independent as possible from the details of the implementation.
We consider the production of a generic colour-singlet system F of invariant mass Q and transverse momentum p T in hadronic collisions. We start with the MiNLO formula [7,9] for an arbitrary infrared-safe observable O, embedded in the POWHEG method [9,28] as follows Equation (2.1) is accurate up to NLO both in the zero and one jet configurations. B denotes the differential cross section for the production of F plus one light parton (FJ), and V and R are the UV-renormalised virtual and real corrections to this process, respectively. The V (Φ FJ ) term in eq. (2.2) is infrared divergent, and so is the integral of R. These divergences cancel in their sum, so thatB is infrared finite. ∆ pwg denotes the usual POWHEG Sudakov form factor, Λ pwg is an infrared cutoff of the order of a typical hadronic scale, and p T,rad corresponds to the transverse momentum of the secondary emission associated with the radiation variables Φ rad .S(p T ) stands for the MiNLO Sudakov form factor [7], that is evaluated using the kinematics Φ FJ in the Born and virtual terms, and with the full real kinematics Φ FJJ in the real term. The factor (1 + αs(pT) 2π [S(p T )] (1) ) is the first order expansion of the inverse of the MiNLO Sudakov form factor, necessary to avoid any source of double counting in eq. (2.1). The MiNLO procedure specifies that the scale at which the strong coupling constant and the parton densities are evaluated should be equal to that contained in the Sudakov form factor, that we take to be the transverse momentum p T of the colour-singlet system F. 4 It also specifies that the transverse momentum appearing in the Sudakov form factor that multiplies R in eq. (2.2) should be the one of the real kinematics configuration Φ FJJ , which differs from the one appearing elsewhere in the formula, that is relative to the underlying Born kinematics Φ FJ . It turns out that in the singular regions of the secondary emission the transverse momentum of F in the real and underlying Born kinematics become identical, so that the cancellation between the collinear and soft singularities can occur.
For simplicity of notation, eqs. (2.1) and (2.2) refer to the case in which there is only one singular region for the secondary emission. In general, there are singularities both in the initial-state (POWHEG handles the two initial-state regions together), and in the final-state. POWHEG deals with the multiple singular regions by partitioning the real matrix elements, as discussed in detail in ref. [28].
To simplify the discussion that follows, without loss of generality, we ignore the first (Sudakov suppressed) term on the right-hand side of eq. (2.1) (this is simply done for the sake of clarity; this term is always included in the POWHEG implementations), and rewrite the equation as In the first two lines of eq. (2.4), the radiation integral has been evaluated according to the usual unitarity condition that is possible because the observable O(Φ FJ ) does not depend upon the radiation phase space Φ rad . Owing to the infrared safety of the observable, the last line of eq. (2.4) has no singularities. This is obvious as far as the secondary emission is concerned, since the difference between the observables vanishes when it becomes unresolved. Singularities associated with the first emission, on the other hand, are suppressed by the fact that the separation of regions in POWHEG [28] ensures that the secondary emission is always more singular than the first one. Therefore, the contribution of the last line of eq. (2.4) is of pure order α 2 s , and it is dominated by large scales. We conclude that in order to achieve NNLO accuracy it is sufficient to correct eq. (2.1) in such a way that it remains unaltered at large p T (where it already has O(α 2 s ) accuracy), and is NNLO accurate for observables of the form O(Φ) = g(Φ F (Φ)), where g is an arbitrary function, and Φ F (Φ) represents an infrared-safe projection of the kinematic configuration corresponding to a generic phase space Φ to the one where the transverse momentum of the colour singlet vanishes (Φ F ). For instance, Φ F involves the rapidity of the colour-singlet system and its internal variables.
According to the MiNLO method, NLO accuracy is guaranteed if theB(Φ FJ ) function in eq. (2.2) is defined as a total derivative up to the relevant perturbative order [7]. As we will show in the next section, achieving NNLO accuracy will require the inclusion of additional terms in theB(Φ FJ ) function.

MiNLO accuracy
In the NNLO+PS approach of ref. [16], NNLO accuracy is achieved by a reweighting procedure. This is carried out by first computing the inclusive cross section at fixed kinematics of the colourless system in the MiNLO approach and at NNLO, and then by reweighting the events by the ratio of the latter result to the former. This procedure works regardless of the corrections that the MiNLO approach already provides at the NNLO level, since this is eventually divided out and replaced by the correct one.
In the present work we are not relying upon a reweighting procedure, and thus we need to develop an analytic understanding of what MiNLO provides at the NNLO order. We do this by noticing that eq. (2.4) is equivalent to the following equation: up to terms of N 3 LO order. This result follows from the fact that the term involving R in the second line of (2.4) cancels exactly the last term in the curly bracket in the last line of (2.4). On the other hand, eq. (2.4) is derived from the full MiNLO result using only the exact unitarity of the shower and of the POWHEG radiation implementation, and thus does not introduce any fixed order approximation. Therefore eq. (2.6) represents analytically what MiNLO provides at the α 2 s level. It has unavoidably a formal character, with the virtual and real contributions that are separately infrared divergent. As such, it is independent of the specific method used to cancel infrared divergences. In particular, it does not depend upon the details of the POWHEG implementation, such as the mapping between the real cross section and the underlying Born, and it can thus be used to make direct contact with analytic resummation formulae.

Reaching NNLO accuracy: the MiNNLO PS method
In this section we present a simple derivation of the missing terms needed to reach NNLO accuracy in the MiNLO formula. For the interested reader, we report a detailed and more rigorous derivation in section 4.
As it will be shown in section 4 (and also appendix E), up to the second perturbative order, the differential cross section in p T and in the Born phase space Φ F is described by the following formula where R f contains terms that are non-singular in the small p T limit. We notice that Φ F on the left-hand side of eq. (2.7) is defined through a projection of the full phase space with multiple emissions, in particular Φ FJ and Φ FJJ , onto the Φ F phase space. We denote this projection by Φ F,res (Φ) , (2.8) and Φ stands for Φ FJ , Φ FJJ , and so on. The suffix "res" in Φ F,res stands for "resummation", to make clear that the projection is relative to how the recoil of the colour-singlet system is treated in the resummation approach. The Sudakov form factorS reads where all coefficients are defined in section 4 and appendix B. The factor L, defined in eq. (4.31) of section 4, involves the parton luminosities, the Born squared amplitude for the production of the colour-singlet system F, the hard-virtual corrections up to two loops and the collinear coefficient functions up to second order. These constitute some of the ingredients necessary for the next-to-next-to-next-to-leading logarithm (N 3 LL) resummation.
Here, for ease of notation, we do not indicate explicitly the Φ F dependence of L and R f .
As it stands, eq. (2.7) is such that its integral over p T between an infrared cutoff Λ (more precisely, the scale value when the Sudakov form factorS(p T ) vanishes) and Q reproduces the NNLO total cross section for the production of the colour-singlet system. We can recast eq. (2.7) as and We now make contact with the MiNLO procedure. We start by writing the regular terms R f to second order as , (2.14) where the notation [X] (i) stands for the coefficient of the i-th term in the perturbative expansion of the quantity X. The first term on the right-hand side of the above equation is the NLO differential cross section for the production of the singlet F in association with one jet J, namely .
(2.15) As a second step, we factor out the Sudakov exponential in eq. (2.11) and obtain (2. 16) We notice that in order to preserve the perturbative accuracy of the integral of eq. (2.16), it is sufficient to expand the curly bracket in powers of α s (p T ) up to a certain order. In fact, when expanded in powers of α s (p T ), all terms in the curly brackets of eq. (2.16) contain at most a 1/p T singularity and (for the terms arising from the derivative ofS) a single logarithm of p T . The contribution of the terms of order α m s (p T ) ln n Q pT to the total integral of eq. (2.16) between the infrared scale Λ and Q is of order [7] Q Λ dp This crucially implies that, for the integral to be NLO accurate, i.e. O(α s (Q)), one has to include all terms up to order α 2 s (p T ) in the curly brackets of eq. (2.16). This guarantees that the perturbative left-over is of formal order O(α 2 s (Q)) in the total cross section. After performing this expansion in eq. (2.16) we obtain . (2.18) We notice that this formula is what we would obtain by integrating formula (2.6) for an observable of the type where Φ F,res (Φ) was introduced in eq. (2.8). It is thus equivalent to the MiNLO formula. This equivalence relies upon the fact that the MiNLO result can be cast in the form of eq. (2.6), that is independent of any particular phase space projection. We stress again that, in order for eq. (2.18) to have NLO accuracy, S must include correctly terms of order up to α 2 s which exactly reproduce the singular part of the cross section and hence ensure that eq. (2.18) can be reassembled back as a total derivative to the desired perturbative order.
In order to achieve NNLO accuracy, it is now sufficient to guarantee that eq. (2.18) has O(α 2 s (Q)) accuracy at fixed Φ F after integration over p T . This requires the inclusion of all terms up to O(α 3 s (p T )) in the curly brackets of eq. (2.16), and we obtain dσ The regular terms that we omitted in eq. (2.20) arise from the O(α 3 s (p T )) expansion of the term R f (p T )/ exp[−S(p T )] in eq. (2.16), which vanish in the limit p T → 0. The absence of a 1/p T singularity ensures that such terms give a N 3 LO contribution to the total cross section, and therefore can be ignored. We explicitly verified that their inclusion yields a subleading numerical effect. Equation (2.20) constitutes the reference formula to build the MiNNLO PS generator. This simply amounts to adding to the MiNLO formula the new term , where all coefficients are defined in appendices B and C.

Implementation of the [D(p T )] (3) term in the MiNLO framework
The MiNLO method based on eq. (2.18) has been implemented within the POWHEG-BOX framework [29] and it has been thoroughly tested. In order to achieve NNLO accuracy, we therefore include the new terms discussed in the previous section as a correction to the existing implementation.
We recall that all terms in the MiNLO formula (2.18) are directly related to the phase space of the production of the colour singlet F together with either one (Φ FJ ) or two jets (Φ FJJ ). Conversely, in the MiNNLO PS master formula (2.20), the new term [D(p T )] (3) arises from a resummed calculation in the p T → 0 limit where the information about the rapidity of the radiation has been integrated out inclusively. As such it depends on the phase space Φ F of the colour singlet with no additional radiation, and carries an explicit dependence on the p T of the system. This dependence, however, does not correspond to a well-defined phase-space point for the full event kinematics (neither Φ FJ nor Φ FJJ ), since the presence of a p T requires at least one parton recoiling against F, but we have no information on the kinematics of such a parton. This has no consequence on the accuracy of the MiNNLO PS formula, since at finite transverse momentum [D(p T )] (3) contributes with a O(α 3 s (Q)) correction to the integrated cross section. It follows that at large values of p T the kinematics associated with the [D(p T )] (3) terms can be completed in an arbitrary way, implying variations beyond NNLO accuracy. In particular, we observe that a Φ F phasespace point can be obtained from a Φ FJ phase-space point through a suitable mapping, while the p T corresponds to that of the Φ FJ kinematics. The mapping should project Φ FJ to Φ F smoothly when p T → 0.
In order to embed the new MiNNLO PS formulation of eq. (2.20) into the MiNLO framework, one must therefore associate each value of [D(p T )] (3) to a specific point in the Φ FJ phase space. This requires supplementing the Φ F and p T information of the [D(p T )] (3) term with the remaining kinematics of the radiation that has been previously lost. In other words, [D(p T )] (3) should be spread over the Φ FJ phase space in such a way that, upon integration, eq. (2.20) is eventually reproduced.
The most obvious way to spread the [D(p T )] (3) term in the Φ FJ phase space is either uniformly, or according to some distribution of choice. To this end, we multiply [D(p T )] (3) by the following factor phase space into the phase space for the production of the colour singlet F alone (for instance performed according to the FKS mapping for initial-state radiation (ISR) discussed in section 5.5.1 of ref. [28], that preserves the rapidity of the colour-singlet system). J is an arbitrary function of Φ FJ , and labels the flavour structure of the FJ production process. Finally, p T is the transverse momentum of the radiation (hence that of the colour singlet) in the Φ FJ phase space. The factor F corr is such that upon integration over the Φ FJ phase space together with a function that depends only on p T and Φ F , the result reduces to the integral of that function over Φ F and p T . In formulae, for an arbitrary function G(Φ F , p T ), we have The full phase-space parametrisation and the Φ FJ → Φ F mapping are given in appendix A. The function J (Φ FJ ) in eq. (3.1) gives us some freedom in choosing how to spread [D(p T )] (3) in the radiation phase space. Among the sensible choices, one could simply use a uniform distribution by setting For this trivial choice an analytic solution for the integral in the denominator of eq. (3.1) is given for illustration in appendix A. However, we found that this choice generates a spurious behaviour when the jet is produced at very large rapidities. A more natural choice is to spread [D(p T )] (3) according to the actual rapidity distribution of the radiation, by setting where |M FJ (Φ FJ )| 2 is the tree-level matrix element squared for the FJ process, and the ) represents the product of the parton densities in the initial-state flavour configuration given by the index . This choice provides a more physical distribution of [D(p T )] (3) , but it can become computationally expensive for complex processes with several degrees of freedom as the integral in the denominator of eq. (3.1) has to be evaluated for every phase-space point numerically. A convenient compromise is to take the collinear limit of the squared amplitude of eq. (3.4), namely where |M F (Φ F )| 2 is the Born matrix element squared for the production of the colour singlet F, and P (Φ rad ) is the collinear splitting function. After noticing that the Born squared amplitude |M F (Φ F )| 2 cancels in the ratio of eq. (3.1), we can simply set where the full expression is reported in eqs. (A.14), (A.13). This prescription is computationally faster, since the integral in the denominator of eq. (3.1) has a better convergence, and it does not change for more involved processes. With these considerations, eq. (2.20) can be recast in a way that is differential in the entire Φ FJ phase space as where the sum over flavour configurations is understood, and p T is meant to be defined in the Φ FJ phase space. A second aspect relevant to the implementation of the MiNNLO PS procedure is related to how one switches off the Sudakov form factor, as well as the terms [S(p T )] (1) and , in the large p T region of the spectrum. We stress that the details of this operation do not modify the accuracy of the result. This is because in the large p T region eq. (3.7) differs from the NLO FJ distribution only by O(α 3 s ) corrections relative to the Born. This implies that one has some freedom in choosing how to turn off the logarithmic terms at scales p T Q. One important constraint to keep in mind is that in the regime p T Q the logarithmic structure has to be preserved in order to retain the NNLO accuracy in the total (inclusive) cross section.
There are of course different sensible ways to switch off the logarithmic terms at large p T . One possibility is to set the quantitiesS(p T ), [S(p T )] (1) , and [D(p T )] (3) to zero at p T ≥ Q. This prescription is adopted in the original MiNLO implementation of ref. [7]. A second possibility, closer in spirit to what is done in resummed calculations, is to modify the logarithms contained inS(p T ), [S(p T )] (1) , and [D(p T )] (3) , so that they vanish in the large p T limit. This is done by means of the following replacement where p is a free positive parameter. Larger values of p correspond to logarithms that tend to zero at a faster rate at large p T . With this modification we also include the Jacobian factor (Q/p T ) p 1 + (Q/p T ) p . (3.9) in front of the [D(p T )] (3) term, which has the effect of switching it off its 1/p T terms when p T goes above Q.
It is easy to convince ourselves that this modification does not alter the MiNNLO PS accuracy. In fact, the modified logarithms only play a role when p T Q. In this regime, the counting of the orders in the MiNNLO PS formula simplifies, since there are no enhancements due to large logarithms. Under these circumstances, the [D(p T )] (3) term is subleading, and the Sudakov form factor also leads to subleading effects once combined with its first order expansion in formula (3.7). Thus, as long as the modified logarithms are used consistently in bothS and [S(p T )] (1) , only terms beyond the relevant accuracy are generated.

Derivation of the MiNNLO PS master formula
In this section we present a more rigorous derivation of the MiNNLO PS formalism that has been outlined in section 2.3. Our starting point is the calculation of the cumulative transverse-momentum spectrum for a colour singlet produced in the collision of two hadrons. More precisely, we consider the second-order perturbative expansion of the above cumulative cross section in the limit p T → 0 (i.e. the singular part), with up to two emissions. This information can be directly accessed in the momentum space formulation of transverse-momentum resummation presented in refs. [26,27]. 5 The singular part expressed in this way reads 6 where we have defined The various terms of eq. (4.2) are explained in the following. We defined the following notation in terms of the initial-state legs a and b: where the identity matrix indicates a trivial dependence on the momentum fraction z, i.e.
Our starting point follows from eqs. (2.58) and (2.59) of ref. [27], by considering only the first two emissions. 6 The convolution between two functions f (z) and g(z) is defined as The regularised splitting functionP [a/b] and the coefficient functions The L factors contain the parton luminosities convoluted with the coefficient functions and multiplied by the virtual corrections. They read (4.11) The last term in eq. (4.2) accounts for azimuthal correlations and it is non-zero only for processes that are gg-initiated at the Born level. Its structure is identical to the one of the first term, provided one replaces the coefficient functions C with the corresponding G functions [30]. The quantity ∆ (C) (Q 1 , Q 2 ) represents the no emission probability between the scales Q 1 and Q 2 < Q 1 , and it is given by and an analogous definition holds for ∆ (G) . All considerations beyond this point hold identically for both the C and G terms, and therefore we omit the latter in the following equations for the sake of simplicity (but its contribution is understood). The function H contains the contribution of the virtual corrections to the colour-singlet process under consideration Note that, when working directly in momentum space, the term H (2) differs from the second-order coefficient of the form factor in the MS scheme [27] (cf. eq. (B.14) below). The various coefficients used in the above equations are reported in appendix B. Finally, the Sudakov radiator S is defined as in eq. (2.9), but with the anomalous dimensionB(α s ) replaced by The first order coefficient B (1) is identical to the one of eq. (2.9), while the difference between the second order coefficients B (2) andB (2) will be explained shortly in this section. All explicit formulae are reported in appendix B. Equation (4.2) reproduces the correct logarithmic structure at O(α 2 s ), including the NNLO constant terms. The evolution, in principle, continues with extra emissions down to the infrared cutoff of the theory Λ Q, that is the scale at which the Sudakov form factor vanishes. However, for the time being we are only considering the first two of such emissions, as the additional radiation will be included by the parton shower via a consistent matching at a later stage.
The formulation that led to eq. (4.2) is based on an organisation of the perturbative series that aims at the resummation of the logarithmic terms to all orders. In particular, each evolution step in eq. (4.2) describes the emission of an inclusive correlated cluster of emissions [27]. In an illustrative picture, it is convenient to think of each inclusive cluster as describing the emission of an initial-state radiation and its subsequent branchings. The cluster is defined by integrating inclusively over the branching variables and retaining only the information about the total transverse momentum of the radiation emitted within. An illustrative example of the separation into correlated clusters for two typical configurations is reported in fig. 1.
In a slightly more technical picture, limiting ourselves to the O(α 2 s ) case we are interested in, the correlated clusters originate from the fact that, in the soft limit, the squared amplitude for the emission of up to 2 partons can be decomposed as The termM 2 (k 1 , k 2 ) is defined as the correlated part of the squared amplitude that cannot be decomposed as the product of two squared amplitudes for the emission of a single parton. 7 We point out that, in the collinear limit, a decomposition of the type (4. 16) requires keeping track of all the possible flavour structures of a 1 → 3 collinear branching. Each of the correlated termsM 2 admits a perturbative expansion in powers of α s defined by including the virtual corrections while keeping the number of emissions fixed. The inclusive correlated cluster is defined as Y ab denotes the rapidity of the k a + k b system in the centre-of-mass frame of the collision, and [dk] denotes the phase-space measure for the emission k. Moreover, the strong coupling constant in eq. (4.17) is evaluated at the transverse momentum of the inclusive cluster. For an inclusive observable such as p T , the only quantity that matters is the total transverse momentum of each inclusive cluster. Therefore, one can integrate over the rapidity of each cluster, and analytically cancel the infrared and collinear singularities. This results in a simplified picture, in which each inclusive cluster of transverse momentum k T,i is emitted according to the kernel and analogously for the term involving the G coefficient function. In eq. (4.18) we implicitly sum over the two initial-state legs a and b, and we observe that the first two terms start at order O(α s ), while the third term proportional to β(α s ), defined in eq. (B.1), starts at O(α 2 s ). Equation (4.2) describes the emission of two subsequent clusters with transverse momenta k T,1 and k T,2 ≤ k T,1 . The first line of eq. (4.2) encodes the emission of the first cluster, while the remaining lines describe the emission of the second one. The latter is split into a no-emission probability term (that excludes configurations with more than one cluster), and a term proportional to the real kernel (4.18) that actually describes the second emission. We stress that, in the present treatment, we are not interested in retaining a given logarithmic accuracy in the p T spectrum, but rather in describing the correct singular structure at O(α 2 s ). However, one should bear in mind that the evolution beyond the second emission will be subsequently generated by a parton shower generator, which will guarantee a fully exclusive treatment of the radiation. We finally observe that the term in the second emission's probability contributes at most to O(α 3 s ) and therefore can be ignored.
We now proceed by adding and subtracting Θ(p T − k T,1 ) to the second Θ function in eq. (4.2) as (4.20) We can then recast eq. (4.2), with O(α 2 s ) accuracy, as where we neglected the Sudakov form factors in the integral containing the difference between the two Θ functions, as the first non-trivial contribution generates O(α 3 s ) corrections. The resulting integral is finite at fixed p T . In order to evaluate such an integral, and since we are only interested in the O(α 2 s ) result, we can expand the content of the curly brackets about k T,1 = p T and k T,2 = p T as follows When performing the integration, we notice that the first term on the right-hand side of the above equation vanishes upon azimuthal integration, and we are left with the integral of the second term, obtaining where S (p T ) = dS (p T )/dL. We now incorporate the terms generated by the latter integration into the first term of eq. (4.23). Retaining O(α 2 s ) accuracy, this can be done by redefining some of the resummation coefficients as follows (an alternative derivation of such redefinitions is performed using an impact-parameter space formulation in appendix E) (4.24) In order to use eq. (4.23) in the context of the MiNLO algorithm, we perform a resummation-scheme transformation to evaluate the virtual corrections (4.14) that appear in L (C) (k T,1 ) (4.11) at a scale k T,1 . This implies the further replacements [7] H(Q) → H(k T,1 ) , (4.25) We thus define and we redefine all ingredients in our calculation as S →S, ∆ (C) →∆ (C) , C →C, and L (C) →L (C) to take the above replacements into account. We therefore recast eq. (4.23) as Finally, we take the derivative in p T in order to obtain the singular structure of the differential p T distribution, that reads In order to be accurate across the whole p T spectrum, we need to match eq. (4.28) to the NLO differential cross section for the production of the colour-singlet system in association with one jet. This can be performed in two steps.
The first step is to observe that the second emission is distributed in a way that closely mimics the treatment of the radiation in the POWHEG method [9] discussed in section 2.1, that is generated according to the probability where the factor Φ rad represents the full FKS [33] radiation phase space for the second emission k 2 . 8 The quantities R and B represent the tree-level squared amplitudes for FJJ (double emission) and FJ (single emission), respectively. Therefore, the second emission can be directly generated according to the POWHEG method, which guarantees an accurate description at tree level for k 2 over the whole radiation phase space Φ rad . We can then focus on the first cluster contribution. For simplicity we can integrate eq. (4.28) explicitly over the second emission k 2 , stressing that the latter can be restored fully differentially by closely following the POWHEG procedure as previously discussed. We obtain where in the second line we recast the result in a more compact form. We can at last restore the contribution of the G coefficient functions by replacingL (C) with the full luminosity factor as (4.31) Equation (4.30), when expanded, correctly reproduces up to O(α 2 s ) the divergent (logarithmic) structure of the differential spectrum in the small p T limit. However, eq. (4. 30) does not yet include the regular terms in the p T distribution (i.e. those which vanish in the p T → 0 limit).
The second step to include the regular terms (i.e. that vanish in the p T → 0 limit) in the above formula is to add the full NLO result for the production of the colour singlet F and one additional jet, and subtract the NLO expansion of the total derivative in eq. (4.30), which leads to where we used eq. (2.14), namely . (4.33) The first term on the right-hand side of eq. (4.32) constitutes the starting point (2.7) for the discussion in section 2.3. In particular, following the considerations that led from eq. (2.7) to eq. (3.7), and restoring the generation of the second radiation via the POWHEG mechanism we obtain where p T is defined in the Φ FJ phase space. Equation (4.34) constitutes the master formula for the MiNNLO PS method, to match a fully differential NNLO calculation to a parton shower. The NNLO subtraction in eq. (4.34) is accomplished thanks to the Sudakov form factor that exponentially suppresses the p T → 0 limit. We stress that this fact does not imply that the transverse-momentum spectrum of the colour singlet will be exponentially suppressed at small p T . The extra emissions beyond the second one, generated by the parton shower, will eventually modify the scaling of the transverse-momentum distribution and restore the correct O(p T ) scaling in this regime [27,34]. This, for sufficiently accurate parton showers ordered in transverse momentum, effectively corresponds to leading logarithmic accuracy in the p T spectrum. 9 5 Application to Higgs-boson and Drell-Yan production at the LHC In this section we apply the MiNNLO PS method to hadronic Higgs-boson production through gluon fusion in the approximation of an infinitely heavy top quark (pp → H), and to the Drell-Yan (DY) process (pp → Z → + − ) for an on-shell Z boson. Rather than presenting an extensive phenomenological study for these two processes at the LHC, our goal is to perform a thorough validation and numerically demonstrate the NNLO+PS accuracy of the MiNNLO PS formula. This requires us to verify two aspects of the results: firstly, that NNLO accuracy is reached for Born-level (Φ F ) observables, in particular for the total inclusive cross section and for distributions in the Born phase space, such as the rapidity distribution of the colour-neutral boson or, in case of DY, the leptonic variables. Secondly, it has to be shown that the NLO accuracy of one-jet (Φ FJ ) observables is preserved, in particular for distributions related to the leading jet. To this end, we compare our MiNNLO PS results to MiNLO and to NNLO predictions. We stress that the results produced by the MiNNLO PS method differ from the previous NNLOPS [16,18] by higherorder terms. In particular, the latter agree by construction with the full NNLO for Born variables, therefore we validate our results by directly comparing to the nominal NNLO. After defining the general setup, we discuss the validation in the following.

Setup
We consider 13 TeV LHC collisions. For the EW parameters we employ the G µ scheme with real Z and W masses, since we consider on-shell Z bosons. Thus, the EW mixing angles are given by The following values are used as input parameters: G F = 1.16639 × 10 −5 GeV −2 , m W = 80.385 GeV, m Z = 91.1876 GeV, and m H = 125 GeV. We obtain a branching fraction of BR(Z → + − ) = 0.0336310 from these inputs for the Z-boson decay into massless leptons. With an on-shell top-quark mass of m t = 172.5 GeV and n f = 5 massless quark flavours, we use the corresponding NNLO PDF set with α s (m Z ) = 0.118 of PDF4LHC15 [36] for Higgs-boson production and of NNPDF3.0 [37] for the DY results.
For MiNLO and MiNNLO PS , the factorisation scale (µ F ) and renormalisation scale (µ R ) are determined by the underlying formalism to be proportional to the transverse momentum of the colour singlet, as discussed in section 2. Upon integration over radiation this corresponds to effective scales (µ R , µ F ) of the order of m H and m Z for Higgs-boson and DY production, respectively. The latter scales are used to obtain the fixed-order results. We stress that the MiNLO and MiNNLO PS methods adopt a dynamical scale during the phase-space integration. As a consequence, the correspondence between such scales and those used in the fixed-order predictions presented below is only approximate, and for this reason one does not expect a perfect agreement between the two calculations.
Uncertainties from missing higher-order contributions are estimated from customary 7-point variations, i.e. through changing the scales by a factor of two around their central values for the fixed-order results) while requiring 0.5 ≤ K F /K R ≤ 2. This implies taking the minimum and maximum values of the cross section for variations (K F , K R ) = (2, 2), (2, 1), (1, 2), (1, 1), (1, 1 2 ), ( 1 2 , 1), ( 1 2 , 1 2 ). The formulae for the scale variation are reported in appendix D. When the scales µ R (µ F ) are too small, we freeze them consistently everywhere, both as arguments of the strong coupling constant and partons densities, and in the terms of the cross section that depend explicitly on them. Our choice for the freezing scale is 1.8 GeV for Drell-Yan, and 2.5 GeV for Higgs-boson production. The choice of these freezing scales is determined from two opposite requirements: one is to stay above the lower bound of the PDF parametrization (which is a scale of the order of 1 GeV), and the other is to remain close to the region where the Sudakov form factor is negligible, in order to make sure that the contribution at the lower integration bound of the quantity defined in (2.7) stays indeed negligible. This allows for a larger value in the Higgs case, where the Sudakov suppression is stronger. We note that our implementation preserves scale compensation and that the transverse momentum is not affected by this procedure. We stress again that the aforementioned scales are only used to freeze the values of µ R and µ F , and they don't act as a cutoff on the phase space Φ FJ .
We notice that, at the freezing scale, the Sudakov form factor is already quite small: in the Higgs case, already at p T = 3 GeV its value is 0.01. In the Drell-Yan case, when p T = 1.8 GeV the value of the Sudakov form factor is 0.1, when p T = 1 GeV it is 0.03, and it reaches the value 0.01 for p T = 500 MeV. We remark that, if we were to exclude from our calculation scale variation points with K F = 1/2, the freezing scale could be taken as low as 1 GeV, which is the PDF cutoff, without a relevant change in the result, as shown in table 1. Finally, we stress that the region affected by the choice of the freezing scale is at transverse momenta where non-perturbative effects start to play a role, and in a realistic simulation these effects are taken into account by the parton shower Monte Carlo and its hadronization model.
It is important to bear in mind that we implement the scale variation in all terms of eq. (4.34), including the SudakovS. The latter variation is not present in a standard NNLO calculation, and therefore probes additional sources of higher-order corrections. Hence, we expect the resulting scale dependence to be moderately larger than the one of the NNLO fixed-order predictions, reflecting the additional sources of perturbative uncertainties in the MiNNLO PS matching procedure. Although such extra sources could be avoided, we reckon that their inclusion is more appropriate to reflect the actual uncertainties of the MiNNLO PS method.
We have employed the POWHEG-BOX framework [9,28,29] [47][48][49][50][51] also with the HNNLO [52] and DYNNLO [53] codes, which we found to be fully compatible with the Matrix predictions within their respective systematic uncertainties. For lepton-related observables we use the results from DYNNLO in our comparison. Unless otherwise stated, all results presented throughout this section are subject to no cuts in the phase space of the final-state particles. All showered results are obtained through matching to the Pythia8 parton shower [54], and they are shown at parton level, without hadronization or underlying-event effects. Finally, the value of the parameter p of the modified logarithms of eq. (3.8) has to be chosen such that the logarithmic terms are switched off sufficiently quickly at large transverse momentum, in order to avoid spurious effects in the region dominated by hard radiation. We adopt p = 6 in the following, that is slightly larger than the values used in standard resummations [27,55], and verified that variations of p lead to very moderate effects that are well within the quoted uncertainties. We report MiNNLO PS results for the total inclusive cross section in table 2, together with the LO, NLO, NNLO, and MiNLO predictions. We start by discussing the Higgs cross sections: Compared to MiNLO we find a +19% effect by including NNLO corrections through the MiNNLO PS procedure. Furthermore, as expected from including an additional order in the perturbative series, the size of the uncertainties due to scale variations are reduced. In particular the upper variation bound is almost a factor of three smaller for MiNNLO PS in comparison to MiNLO . The MiNNLO PS and NNLO results agree well within their respective scale-uncertainty bands, which largely overlap. The central values of the two calculations differ by 7.9% as can be seen from the ratio, and they are included in the uncertainty bands of each of the two calculations. The size of the difference is justified by the large perturbative corrections that characterise Higgs-boson production, which implies that subleading terms can be sizable. For processes with smaller corrections these differences will reduce, as we will see in the context of the DY process below. As already mentioned above, it is not expected that MiNNLO PS reproduces the NNLO result exactly, as subleading corrections are treated differently in the two calculations and the the renormalisation and factorisation scales are set differently. Consistency between the two NNLO-accurate predictions within perturbative uncertainties is therefore sufficient, and shows that the MiNNLO PS procedure induces the expected corrections. The scale uncertainties of MiNNLO PS for the total cross section (∼ 13%) are slightly larger than the NNLO ones (∼ 10%). There are two main reasons for this behaviour: On the one hand, MiNNLO PS probes scales in both the PDFs and in α s that in the bulk-region of the cross section are much lower (∼ p T ) than in the fixed-order computation, which naturally induces a larger scale dependence. On the other hand, we include additional scale-dependent terms (as pointed out before) that originate from the analytic Sudakov form factor in the MiNNLO PS procedure, which are absent in a fixed-order calculation, see appendix D. This induces a more conservative estimate for the theory uncertainties of the MiNNLO PS predictions.

Inclusive cross section
In the case of the DY results in table 2, we observe that conclusions similar to the case of Higgs production can be drawn, albeit with significantly smaller corrections: The effect of the MiNNLO PS procedure is to increase the MiNLO cross section by about 5%. Again the scale uncertainties are vastly reduced, in the case of DY by almost a factor of 10. The MiNNLO PS result is only 1.7% below the NNLO prediction and they are in good agreement within their respective scale uncertainties, which are extremely small. Roughly speaking, scale uncertainties are 2% for MiNNLO PS , which is a bit larger than the 1% uncertainties at NNLO. Given the above discussion about the formal differences between MiNNLO PS and NNLO fixed-order computations, these results are very compelling and provide a numerical proof of the accuracy of the total inclusive cross section of the MiNNLO PS procedure. We will now turn to validating the MiNNLO PS results also for differential observables.

Distributions for Higgs-boson production
We first consider the case of Higgs-boson production. The figures of this section are organized as follows: the main frame shows the results from MiNNLO PS (blue,solid) and MiNLO (black, dotted) after parton showering, as well as NNLO predictions (red, dashed), and all results are reported in units of cross section per bin (namely, the sum of the values of the bins is equal to the total cross section, possibly within cuts). In an inset we display the bin-by-bin ratio of all the histograms which appear in the main frame to the MiNNLO PS curve. The bands correspond to the residual uncertainties that are computed from scale variations as indicated in section 5.1.
The transverse-momentum distribution of the Higgs boson (p T,H ) is shown in the left panel of figure 2. At fixed order this distribution diverges in the p T,H → 0 limit, and the accuracy is effectively reduced to NLO across the spectrum. By comparing MiNNLO PS and MiNLO curves, we observe that the NNLO corrections are included consistently in  figure 2 is the most relevant observable for which MiNNLO PS needs to be validated against the NNLO result. Indeed, we find that up to statistical fluctuations the NNLO/MiNNLO PS ratio of the distribution is completely flat, which shows their equivalence. Henceforth, the difference of the two results is purely due to the normalisation, i.e. the total inclusive cross section, which has been discussed in detail in section 5.2 and requires no further comments. In particular, the conclusions about the uncertainty bands and the size of the corrections drawn from table 2 hold also for the rapidity distribution shown in figure 2.
We conclude our discussion of the results for Higgs-boson production by looking at jetrelated distributions. We note that the transverse-momentum distribution of the leading jet is very similar to the one of the Higgs boson which is why we refrain from showing it here and refer to the discussion for p T,H . Figure 3 shows the rapidity distribution of the leading jet in the left panel, and the rapidity difference between the leading jet and the Higgs boson in the right panel. Jets in this case are defined using the anti-k T clustering [56] with a radius R = 0. result does not differ from the MiNLO prediction significantly, i.e. beyond perturbative uncertainties. In particular, this numerical check is important to ensure that the implementation and spreading of the [D(p T )] (3) terms in the Φ FJ phase space, as described in section 3, is appropriate. We refrain from showing the NNLO curve for these distributions, as it does not add any relevant information to these tests. Indeed, we observe that, by and large, the MiNLO /MiNNLO PS ratio is flat for both distributions in figure 3, and that the two results agree very well within perturbative uncertainties. We have repeated these checks for various p T,J thresholds in the jet definition, with the same conclusions. In particular, we found that for hard configurations (p T,J 60 GeV) the MiNNLO PS and MiNLO results become essentially identical, as expected from the fact that MiNNLO PS induces no additional corrections in phase-space regions where the radiation is hard. Furthermore, a similar level of agreement is found also for the azimuthal angle between the leading jet and the Higgs boson.

Distributions for Drell-Yan production
We now move on to discuss distributions for the DY process. Since most of the conclusions are similar to the ones for Higgs-boson production, we keep the discussion rather brief. In addition to the results discussed for the Higgs, we also study the kinematics of the leptons arising from the decay of the Z boson. Figure 4 shows the transverse-momentum distribution of the Z boson (p T,Z ) in the left panel, and its rapidity distribution (y Z ) in the right panel. As seen before, the corrections are smaller in the case of the DY process, but the general behaviour is the same as for Higgs-boson production: At large p T,Z the MiNNLO PS result is essentially identical to the MiNLO one, while the additional NNLO terms enter at smaller values of p T,Z . The  NNLO spectrum diverges at small p T,Z , and is harder in the tail due to the different scale setting. For the y Z distribution, the MiNNLO PS uncertainties are significantly reduced with respect to the MiNLO ones. In the central region (|y Z | < 3) the NNLO/MiNNLO PS ratio is nicely flat up to statistical fluctuations, and the two results agree within their respective uncertainties. For very forward Z bosons (|y Z | > 3), on the other hand, we observe a slight increase of the NNLO/MiNNLO PS ratio. We have checked explicitly that without the Pythia8 parton shower, i.e. at the level of Les Houches events, this effect is more moderate and the NNLO and MiNNLO PS uncertainty bands overlap in the forward region. In fact, we noticed that already for the MiNLO prediction, Pythia8 has the same effect, making the Z-boson rapidity distribution slightly more central. 10 Next, we consider the transverse-momentum distribution of the negatively charged lepton (p T, − ) and its rapidity distribution (y − ) in the left and right panels of figure 5, respectively. For the rapidity distribution the relative behaviour between MiNNLO PS , MiNLO , and NNLO is essentially identical to the one of the Z-boson rapidity and does not require any further discussion. As far as the transverse-momentum spectrum is concerned, the NNLO result shows a very peculiar behaviour for p T, − = m Z /2, which reflects the perturbative instability associated with the fact that the leptons at LO are back-to-back and can share only the available partonic centre-of-mass energy √ŝ = m Z , so that their transverse momenta can be at most p T, − ≤ m Z /2. Beyond this value the NNLO result is therefore effectively only NLO accurate, which can be also seen from the increased uncertainty band. Since such an instability is related to soft-gluon effects, this feature is cured in both the 10 We observed that part of this effect can be attributed to the global recoil adopted by Pythia8 for ISR.
The difference from the NNLO prediction is reduced if one uses a more local scheme for the parton-shower recoil, e.g. via the flag SpaceShower:dipoleRecoil=1. MiNNLO PS and MiNLO results, which are in good agreement with each other in terms of shape. Again the MiNNLO PS uncertainty band is significantly smaller than the MiNLO one, and we observe a rather constant correction, of the order of ∼ 5 − 10%, due to the additional NNLO terms. Finally, also for the DY process the jet-related observables are fully consistent within uncertainties when comparing MiNNLO PS and MiNLO predictions, as can be seen in figure 6. However, the size of their uncertainty bands is very different. This is due to the fact that in the original MiNLO prediction a different prescription for the scale variation was adopted, that also involved the integration boundaries of the Sudakov form factor. We have checked that by using our prescription in MiNLO the uncertainty band becomes comparable to the MiNNLO PS one. We stress again that we have tested a variety of p T,J thresholds in the jet definition, and also looked at the azimuthal angle between the leading jet and the Z boson, and found consistent results throughout.

Summary
In this article we have presented a novel approach, dubbed MiNNLO PS , to combine NNLO QCD calculations with parton showers for colour-singlet production at the LHC. The method is based on the MiNLO procedure, which achieves NLO accurate predictions simultaneously in the zero-jet phase space Φ F and in the one-jet phase space Φ FJ . The necessary terms to achieve NNLO accuracy are derived by establishing a connection of the MiNLO and POWHEG methods with the structure of transverse-momentum resummation in direct space. The consistent inclusion of these terms on top of a MiNLO computation allows us to achieve NNLO+PS accuracy for a variety of collider reactions.
We have discussed in detail a suitable implementation of the NNLO corrections within the MiNLO formalism, and their spreading in the Φ FJ phase space. The resulting matching preserves the leading logarithmic structure of the shower Monte Carlo for showers ordered in the transverse momentum, and the final result is NNLO accurate in the zero-jet phase space while being NLO accurate in the one-jet phase space. The combination of the two multiplicities does not require any unphysical merging scale.
As a proof of concept, we have applied the approach to hadronic Higgs production in the heavy-top limit and to the DY process, where a pair of leptons is produced via the decay of an on-shell Z boson. Our results show that NNLO accuracy is reached both for the total inclusive cross section and for Born-level distributions. Differences with NNLO fixed-order results arise only from terms beyond the nominal accuracy, and the two calculations agree well for such observables within the respective perturbative uncertainties estimated from scale variations. As expected, we observe a significant reduction of the scale dependence with respect to the MiNLO results, in line with the inclusion of the NNLO corrections. It was further verified that for jet-related observables in the Φ FJ phase space, where the accuracy of MiNNLO PS and MiNLO is formally identical, no significant effects are induced by the MiNNLO PS corrections.
The algorithm is very efficient, and NNLO accuracy is achieved directly at generation time without any additional reweighting. The total MiNNLO PS simulation requires just 50% more CPU time than the usual MiNLO computation. This makes it suitable for the application to more involved colour-singlet processes, such as vector-boson pair production, which is of significant phenomenological interest. A potential limitation of the algorithm concerns systems with very low invariant mass, such as low-mass diphoton production, whose p T distribution is peaked towards non-perturbative scales. In this situation, the Sudakov form factor, which is responsible for the NNLO subtraction of the infrared singularities, can become intrinsically non-perturbative. Nevertheless, such scenarios are commonly not of experimental interest. Finally, the MiNNLO PS approach could be generalized to the production of massive coloured final states, such as top-quark pair production. Detailed studies of further applications of the MiNNLO PS method are left for future work.

A Phase-space parametrisation for the [D(p T )] (3) term
In this appendix we define the phase-space mapping from Φ FJ to Φ F adopted for initialstate radiation in POWHEG, and discussed in section 5.5.1 of ref. [28]. The projection is defined by performing a longitudinal boost of the FJ system to a frame where F has zero rapidity, followed by a perpendicular boost that modifies the transverse momentum of F so that it is equal to zero, followed by a longitudinal boost, exactly opposite to the first one, that restores the original rapidity of F. After this sequence of boosts, the rapidity of F remains unchanged, but its transverse momentum has become zero, thus yielding a kinematic configuration in the Born phase space Φ F . The Φ FJ phase space can then be expressed in a factorised form: where s is the square of the total incoming energy, and which are all defined in the centre-of-mass frame of the FJ system. The transverse momentum is given by So, we get (A.6) We now replaces = s(1 − ξ), wheres is the virtuality of the FJ system, and multiply and divide by 2p T . After rearranging the delta function we get which is a monotonically increasing function of ξ in the range [0, ∞] for ξ ∈ [0, 1]. We have and obtain and t max is defined as the t value corresponding to ξ max , namely , (y → −y, 1 → 2) .
The variablesx 1,2 are the momentum fractions of the two initial-state partons in the Φ F phase space, and are defined as [28] x with x 1,2 being the momentum fractions of the two initial-state partons in the Φ FJ phase space. For most choices of the function J (Φ FJ ) discussed in section 3, the above integral must be evaluated numerically via usual Monte Carlo techniques. However, for the choice J (Φ FJ ) = 1, we can perform the y integration analytically. We need the following elementary integral . (A.11) The limits of integration have to be computed numerically, by finding the minimum and maximum value such that the Θ function in eq. (A.9) is 1. Denoting by y 1 and y 2 the lower and upper integration boundaries, respectively, we simply find that Finally, we conclude this appendix by reporting the collinear approximation for J (Φ FJ ) discussed in section 3. Its expressions are taken from ref. [57] and adapted to the notation of this appendix. For gluon-initiated processes, the different flavour configurations read For quark-initiated reactions we have (A.14)

B Resummation formulae
In this section we report the expressions of the quantities appearing in the calculation of the analytic transverse-momentum spectrum that we have used throughout this article. First of all we report our convention for the renormalisation-group equation of the strong coupling: where the coefficients of the β-function are 2Nc , N c = 3, and the number of light flavours n f = 5. The Sudakov radiatorS(p T ) in eq. (2.9), with the accuracy considered in this article, can be expressed asS with λ = α s (Q)β 0 ln(Q/p T ), and The resummation coefficientB (2) is defined as according to eqs. (4.24), (4.25), namelỹ For Higgs-boson production in gluon fusion, the coefficients A (i) and B (i) which enter the formulae above are Similarly, for Drell-Yan production they read The expressions for the coefficients A (i) and B (i) are extracted from ref. [58,59] for Higgsboson production and ref. [60] for DY production. The hard-virtual coefficient functions H andH up to two loops are given by for the DY process. The extra term is a feature of performing the resummation in momentum space, and does not appear in the impact-parameter (b) space formulation of transverse-momentum resummation (see ref. [27] for details). The coefficientH (2) , that appears in eqs. (4.26) and (4.31) reads Finally, we report the expansion of the collinear coefficient functions C ab ,C ab , G ab where µ is the same scale that enters parton densities. The first-order expansion has been known for a long time and reads The second-order collinear coefficient functions C (2) ab (z), as well as the G coefficients for gluon-fusion processes are obtained in refs. [61][62][63], while for quark-induced processes they are derived in ref. [64]. In the present work we extract their expressions using the results of refs. [61,64]. For gluon-fusion processes, the C (2) gq and C (2) gg coefficients normalised as in eq. (B.16) are extracted from eqs. (30) and (32) of ref. [61], respectively, where we use the hard coefficients of eqs. (B.12) without the momentum-space term ∆H (2) in the expression for the H (2) (Q) coefficient. 11 The coefficient G (1) is taken from eq. (13) of ref. [61]. Similarly, for quark-initiated processes, we extract C (2) qg and C (2) qq from eqs. (32) and (34) of ref. [64], respectively, where we use the hard coefficients from eqs. (B.13) without the momentum-space term ∆H (2) in the expression for the H (2) (Q) coefficient. The remaining quark coefficient functions C (2) qq , C (2) qq and C (2) qq are extracted from eq. (35) of the same article. The coefficientC (2) (z), that appears in eqs. (4.26) and (4.31) finally reads Accordingly, eq. (D.1) becomes dσ dΦ F dp T = d dp T exp[−S(p T )]L(Φ F , p T , K R , K F ) + R f (Φ F , p T , K R , K F ) , (D. 5) and includes all relevant scale-dependent terms at NNLO. Formula (D.5) retains its NNLO accuracy whether or not we include also scaledependent terms in the Sudakov form factorS. However, in order to make contact with the POWHEG formula, the scale dependence inS must be included. In fact, if we take the derivative in eq. (D.5) we obtain dσ dΦ F dp T = exp[−S(p T )] − dS(p T ) dp T L(Φ F , p T , K R , K F ) + dL(Φ F , p T , K R , K F ) dp T Since in POWHEG we do not have access separately to the terms arising from the derivative ofS, and all terms in the curly bracket are evaluated with the same scale choice, we must make sure that also the derivative ofS is given in terms of α s (K R p T ). This is achieved by writingS asS (p T ) = 2 Q pT dq q A(α s (K R q), K R ) ln Q 2 q 2 +B(α s (K R q), K R ) , (D. 6) where the A and B coefficients include scale-compensating terms in such a way that they are formally independent upon K R when summed up to all orders in perturbation theory. It is easy to see that, with this replacement, the form ofS given in eq. (B.4) remains the same provided that the A (i) and B (i) coefficients are replaced by the K R dependent ones, and that α s (Q) is replaced by α s (K R Q). We now present in detail the formulae needed to implement the scale variation. We start by discussing the L factor, defined in eq. (4.31). The coefficients H (1) andH (2) become H (2) (K R ) =H (2) + 4 n B 1 + n B 2 π 2 β 2 0 ln 2 K 2 R + π 2 β 1 ln K 2 R + 2 H (1) (1 + n B ) πβ 0 ln K 2 R , (D.7) with n B being the α s power of the Born cross section for the production of the colour singlet F. The coefficient functions C receive the following scale dependence: C (1) (z,K F ) = C (1) (z) −P (0) (z) ln K 2 F , C (2) (z,K F , K R ) =C (2) (z) + πβ 0P (0) (z) ln 2 K 2 F − 2 ln K 2 F ln K 2 R −P (1) (z) ln K 2 F + 1 2 (P (0) ⊗P (0) )(z) ln 2 K 2 F − (P (0) ⊗ C (1) )(z) ln K 2 F + 2πβ 0 C (1) (z) ln K 2 R , (D. 8) while G (which is present only in the case of gluon-induced reactions) remains unchanged. We then consider the Sudakov radiatorS, defined in eq. (2.9). We change the scale of the strong coupling in its integrand (2.9) from p T to K R p T , and modify the A and B coefficients as follows 12 The term proportional to n B (the power of α s at LO), is induced by the presence of H (1) in theB (2) coefficient, that in turn originates from evaluating the hard virtual corrections at p T in the factor L, see eq. (4.25).
The scale dependence also propagates into the constituents of the [D(p T )] (3) term, whose α 3 s prefactor in eq. (4.34) is evaluated at K R p T . Besides the dependence in the coefficients reported above (which is understood in the equation that follows), [D(p T )] (3) acquires additional explicit scale-dependent terms:

E Considerations from impact-parameter space formulation
In this section, we derive the form of the starting equation (2.7) using the impact-parameter space formulation of transverse-momentum resummation. We start from the formula and λ b = α s (Q)β 0 ln(Qb/b 0 ), b 0 = 2e −γ E . The g i functions are analogous to those used in momentum space (B.7), and [65] The factor L b is defined as whereH is identical to H of eq. (B.11), with the only difference being that theH (2) coefficient does not contain the term ∆H (2) (B.14). We evaluate the b integral by expanding b 0 /b about p T in the integrand. While this procedure is known to generate a geometric singularity in the p T space resummation, in this article we are only interested in retaining O(α 2 s ) accuracy and therefore this is not an issue for the present discussion. We follow the appendix of ref. [55], and by neglecting terms that contribute beyond O(α 2 s ), we obtain where S (p T ) = dS (p T )/d ln(Qb/b 0 ). After performing the derivatives, we observe that, retaining O(α 2 s ) accuracy, we can approximate the above equation as follows We directly observe that the two terms proportional to S are analogous to those produced in the last line of eq. (4.23). These two terms can be incorporated in the master formula via the replacements (4.26). On the other hand, the term proportional to S is a new feature of the b-space formulation, and it is not present in the momentum space formulation. Using the expression S (p T ) = 32A (1) πβ 0 α 2 s (2π) 2 + O(α 3 s ) (E. 7) we observe that the new O(α 2 s ) constant term 8/3 ζ 3 A (1) πβ 0 can be absorbed into the coefficientH (2) asH (2) → H (2) =H (2) + 8 3 This is precisely the difference between theH (2) coefficient (defined in b space) and the H (2) coefficient present in the momentum-space formulation. Therefore, the O(α 2 s ) expansion of eq. (E.6) coincides with that of eq. (2.7).