Extending the MINLO method

We consider improving POWHEG+MINLO simulations, so as to also render them NLO accurate in the description of observables receiving contributions from events with lower parton multiplicity than present in their underlying NLO calculation. On a conceptual level we follow the strategy of the so-called MINLO' programs. Whereas the existing MINLO' framework requires explicit analytic input from higher order resummation, here we derive an effective numerical approximation to these ingredients, by imposing unitarity. This offers a way of extending the MINLO' method to more complex processes, complementary to the known route which uses explicit computations of high-accuracy resummation inputs. Specifically, we have focused on Higgs-plus-two-jet production (HJJ) and related processes. We also consider how one can cover three units of multiplicity at NLO accuracy, i.e. we consider how the HJJ-MINLO simulation may yield NLO accuracy for inclusive H, HJ, and HJJ quantities. We perform a feasibility study assessing the potential of these ideas.


Introduction
In recent years, next-to-leading order parton shower (Nlops) matching techniques have been developed and realized as practical simulation tools, routinely used in LHC data analysis [1][2][3][4][5][6][7]. By now Nlops methods have been applied to many processes involving the production of a primary colourless system, e.g. a massive boson, B, in association with jets (Bnj) [8][9][10][11][12]. A Bnj Nlops simulation yields NLO accuracy for B + n-jet inclusive observables, and LO accuracy for B + m-jet ones (m = n + 1), while its predictions for more inclusive observables are divergent. Motivated by the success of leading order matrix element-parton shower multi-jet merging approaches in the earlier part of the last decade [13][14][15][16], it has been considered highly desirable to combine Nlops generators for Bnj processes corresponding to different jet multiplicity, n, to obtain a unified simulation output, consistently describing inclusive B, B + 1-jet (Bj) and B + 2-jet (Bjj) observables, simultaneously, with NLO accuracy.
All of these approaches separate the output of each component simulation (B, Bj or Bjj) according to the jet multiplicity of the events it produces, discarding those having a multiplicity for which the generator does not possess the relevant NLO corrections. Having processed the output of each simulation in this way, the event samples are joined to give an inclusive sample. Loosely speaking, each generator can be regarded as contributing a single exclusive jet bin to the final inclusive sample, the magnitude of each bin being predominantly determined by the jet resolution scale used in performing the merging, the so-called merging scale. Different approaches use different means to mitigate the dependence on this unphysical scale.
If the merging scale is too high one loses the benefits of the higher multiplicity generators, describing relatively hard jets with tree-level accuracy, or the parton-shower approximation. If the merging scale is too low, the inclusive sample is dominated by the higher-multiplicity generators, which in general leads to unitarity violation, whereby more inclusive quantities like the total inclusive cross section, exhibit spurious differences with respect to their corresponding conventional NLO predictions. The Geneva approach [22] can completely avoid unitarity violation, and even the introduction of a merging scale, by employing very high accuracy resummations. The method has been demonstrated for effectively merging two units of multiplicity (without a merging scale) in the context of an NNLL +NNLO parton shower matched simulation of the Drell-Yan process [24].
In the sense that it proposes to resolve the merging problem through implementing sufficiently high order resummation, Geneva represents the best solution of merging problem. However, details pertaining to exactly how this is done are subject to debate in the community. In the Unlops approach [25] unitarity is exactly maintained for sufficiently inclusive quantities, through what the authors refer to as 'subtract-what-you-add' approach. Nevertheless, the Unlops method is affected by other complications connected to the presence of a merging scale.
In the Minlo framework [26], fully differential NLO cross sections for processes of type B + n − jets are matched onto a leading log resummation of the exclusive n−jet cross section, as defined by the k t -jet algorithm, with B here referring typically to a given collection of colourless final-state particles. This is the generalization, to the NLO level, of the resummation applied in the Ckkw formalism [14,27] to the highest multiplicity tree level matrix element [16]. Essentially the nhardest pseudopartons found by the k t -jet algorithm have a distribution which is equivalent to that of a parton shower simulation of B-production with, in addition, matching to the exact NLO Bnj matrix elements.
Whereas previously the leading order parts of these cross sections themselves would exhibit unphysical divergences when the Born partons became soft and/or collinear to one another, with the Minlo prescription applied their behavior is instead regular, physical and sensible, i.e. Bnj computations with the Minlo prescription also yield physical results for Bmj (m < n) and even fully inclusive B-production observables. In the case of Bj, with B a W/Z/Higgs boson, it was found that the standard Minlo procedure yielded results for inclusive B production observables equivalent to conventional NLO ones up to terms O(α 3/2 S ) relative to the LO component [28]. In the same article it was proven how, by delicate adjustment of the Minlo Sudakov form factor and clustering procedure, the spurious O(α 3/2 S ) terms could be eliminated, with the subsequent Bj-Minlo calculations achieving NLO accuracy for both B and Bj inclusive observables; in the following we call the improved Minlo procedure of ref. [28] Minlo to distinguish it from the original Minlo prescription [26]. Thus one obtains, from the single NLO calculation of Bj production, also the fully differential NLO calculation of B production. The Minlo calculation can then be matched to a parton shower using the standard techniques [1][2][3]. When viewed in the context of the recent work in merging multiple NLO calculations together this amounts to merging without any unphysical merging scale. It was also demonstrated in refs. [29][30][31] how to promote the Minlo simulations to Nnlops simulations.
While the modifications made in going from Minlo to Minlo involve including higher order terms in the Sudakov form factor, and lead to the recovery of NLO accuracy also for inclusive B observables, the related resummation is not improved in accuracy. The resummation of the B system's transverse momentum is NNLL σ accurate 1 [32,33] before and after the inclusion of the latter terms in the Sudakov form factor [28] (and before Nlops matching adds ambiguities).
Thus, Minlo amounts to the Minlo method with additional, subtle, unitarization. This is in difference to the Geneva approach, wherein higher order resummation is taken as the main defining specification, with unitarization coming 'for free' along with the latter [22,24,34]. In this sense Minlo is, minimally, the same as the Powheg method. Indeed, in the Powheg method a very specific Sudakov form factor is required to achieve an exact unitarization, needing terms in the exponent that are sub-leading with respect to the resummation accuracy to do so, including even power suppressed terms that are nothing to do with resummation.
To realize the Minlo method one needs to know the v → 0 singular part of the Bj cross section differential in the underlying Born variables, Φ B , describing the kinematics of the produced B-final state, where v is a variable measuring radiation hardness, at NLO; this information may be obtained from suitably integrated NLO predictions for the spectrum, or from fixed order expansion of N 3 LL σ resummation. In the case of ref. [28] v was given by the transverse momentum of the produced W/Z/Higgs boson, for which the latter NLO distributions have long been known in the literature [35][36][37]. Recently these distributions have also become available to the same level of accuracy for the transverse momentum of the hardest produced jet [38][39][40][41][42][43][44], with which an equally accurate Minlo calculation could have been developed. For more complicated observables, such as those which might be used for implementing the Minlo method in the context of highermultiplicity processes, with the exception of the N -jettiness variable [45][46][47][48][49][50][51][52], these distributions are (so far) not available in the literature. We note, however, that important progress is being made in the direction of automated approaches to final-state resummation at NNLL [53], valid for broad classes of observables, including those we consider in the present work. Whenever these theoretical ingredients become available, the Minlo method is in principle straightforward to apply, to make a NLO Bnj calculation simultaneously NLO accurate for Bmj (m = n − 1) observables etc; many of the details for that are clarified by the present work.
Nevertheless, even when all the necessary theoretical ingredients are at hand, experience with the Bj-Minlo calculations tells that the results of implementations are sometimes not as ideal as one might have liked. The Minlo codes are proven to return conventional NLO results for inclusive B and Bj observables up to NNLO sized ambiguities and power corrections. In practice the Hj-Minlo calculation was found to give very satisfactory agreement with the regular NLO predictions for inclusive Higgs boson production. On the other hand, comparing Wj-Minlo and Zj-Minlo predictions for inclusive W and Z production to those of regular NLO calculations, one can see numerical differences between the two sets of formally equivalent results, which don't really sit easily with the fact that the two formally agree up to NNLO-sized ambiguities. In the Wj-Minlo case, the inclusion of the relevant NNLL terms in the Sudakov form factor do not lead to noticeably better agreement with the conventional NLO cross sections than those obtained with the original Minlo prescription. We also point out that the true NNLO corrections to Higgs boson production are large, ∼ 20 − 30%, thus the almost perfect agreement of Hj-Minlo with conventional NLO calculations for inclusive Higgs boson production -which looks to be a striking vindication of the theoretical framework -should be considered fortuitous.
Some people (like us) may be dismissive of numerical disagreement between Bj-Minlo and standard NLO predictions for fully inclusive observables, since the Minlo method has been rigorously proven. Others may be less comfortable accepting the fact that these differences arise from contributions beyond the formal accuracy of either type of calculation, given their size in some cases. If 5 − 6% differences can be found in total inclusive cross sections for inclusive W and Z production, it does not seem unreasonable to expect that larger differences may be found in more complex processes, with more powers of the strong coupling associated to the LO contribution and a richer kinematic content. Assuming one is content to dismiss differences due to higher order ambiguities, for complex processes, with even more complex calculations underlying them, it will be difficult to satisfy oneself that the level of numerical disagreement is of this kind.
A final motivation for considering to extend the reach of the Minlo method is that of merging NLO calculations differing by more than one unit of jet multiplicity. Specifically, one would ultimately like a Minlo procedure applied to Bjj-Minlo such that it retrieves NLO accuracy for inclusive B, Bj and Bjj observables. Naive extension of the Minlo method then implies having a N 3 LL σ -accurate nested resummation with which to base it on. While the resummation community is making impressive progress in recent years [54], the prospects for obtaining the high-accuracy ingredients needed to tackle this issue in the near future are unclear to us.
Noting these desirable and undesirable features of the existing Minlo method, we investigated extending it in a number of ways: 1. The Minlo specification can be reached with only limited knowledge of the required Sudakov form factors (at least NLL σ ).
Thus, one can begin to make Minlo simulations for more complex processes.
2. Bjj/Bj-Minlo predictions for Bj/inclusive B observables agree precisely with those of the corresponding conventional NLO calculations.
Numerical ambiguities between conventional inclusive NLO predictions and the associated Minlo ones can be largely eliminated. method can also be used for the purposes of rendering Bjj Minlo simulations simultaneously NLO accurate in the description of inclusive B and Bj production. Section 4 presents a feasibility study assessing the potential of these ideas. We summarize our findings and conclude in sect. 5.

Merging two units of multiplicity
In the following we ultimately present a method for merging Nlops simulations of B-and Bjproduction and, separately, Bj-and Bjj-production without any actual merging. More precisely, the improved Minlo procedure will render the Bj simulation also NLO accurate for B-production and, in the case of Bjj it will build in NLO precision for Bj.
We remind that in this work we refer to the leading tower of logarithms in the spectrum, terms O(ᾱ n S L 2n−1 ), as LL σ , with NLL σ denoting the next-to-leading log tower, O(ᾱ n S L 2n−2 ), NNLL σ for the next-to-next-to-leading log tower O(ᾱ n S L 2n−3 ), and so forth. In section 2.1 we introduce preliminary notation and definitions, in particular regarding the clustering variables which the Minlo procedure is to resum, and the so-called underlying Born kinematics that the resummation is performed about. In section 2.2 we present a formula for NNLL σ resummation of these clustering variables (k t -jet resolutions) based on the Caesar resummation framework [55]. The Sudakov form factors of the latter are compared to the corresponding Minlo formulae in section 2.3. In section 2.4 we derive the fixed order expansion of the NNLL σ Caesar formula and from this we show how our Minlo procedure applied to the Bj(j) NLO computations returns a matched, resummed, NLO accurate jet resolution spectrum. In doing so, we also assume that spurious, unknown, NNLL σ or N 3 LL σ terms may arise in the Minlo resummation, owing to a lack of understanding of the true N 3 LL σ spectrum truncated at NLO, and we closely monitor how these propagate through the Minlo procedure, to better understand and eradicate them, as needed. In section 2.5 we integrate over the Minlo jet resolution spectrum, determining that the distribution of the Born kinematics differs from that of conventional NLO owing to the latter spurious terms, which we have tracked and quarantined. In section 2.6 we first demonstrate how the original Minlo approach removes such terms, by explicitly correcting the Minlo Sudakov form factor, leading to NLO accurate Born kinematics. In the second part of section 2.6, we present our new proposal. Over-simplifying somewhat, this amounts to using the constraint that the corrected Minlo predictions must recover NLO (or NNLO) results for the inclusive Born kinematics, as a defining equation for the elusive Sudakov correction factors, needed to promote Minlo to Minlo .
This equation can no doubt be solved for the latter in several ways, and we present one basic, simple, way to do so. In order to maintain NLO accuracy in the higher multiplicity phase space, it is necessary that the initial Minlo resummation be at least NLL σ accurate, however, this is an easily obtainable threshold by today's standards.
We underline now that it is a working assumption of the Caesar framework, that the underlying Born kinematics, about which the resummation of soft radiation is carried out, are not themselves associated with large logarithmic corrections, i.e. it is assumed that the radiating particles in the hard underlying Born kinematics are well separated. For the case of B production, with only two radiating particles in the initial-state, the latter criterion is fulfilled automatically. On the other hand, for Bj production it implies that we have to restrict ourselves to a regime in which the finalstate (pseudo-)parton in the underlying Born has transverse momentum of order the mass of B, or greater. In other words, in this section 2 it should be understood that the y 12 resummation in Bj-production assumes y 01 O(m 2 B ). Only in section 3 will we consider extending down into the region where the transverse momentum of the final-state Born (pseudo-)parton is small.

Definitions: jet resolutions and underlying Born kinematics
Since it is underlies the whole discussion we first quickly present a reminder of the exclusive k t -jet clustering procedure (for brevity we henceforth refer to pseudopartons obtained in the clustering sequence as just partons): 1. In an n-parton final-state we determine the smallest distance obtained by considering all pairwise combinations of final-state partons i and j, with k t,i , y i and φ i being, respectively, the transverse momentum, rapidity and azimuth of parton i, and {d iB } the set of all final-state transverse momenta: R is the so-called jet radius parameter which we take equal to be 1.
2. If d (n) = d ij partons i and j are replaced in the event by a single mother parton with four momentum p i + p j , otherwise, if d (n) = d i , particle i is considered to have been similarly absorbed in one of the beam jets and is deleted from the event.
3. If further partons remain return to step 1, otherwise the clustering sequence terminates.
In order to have the Bj-Minlo calculation return NLO accuracy for inclusive B production and Bjj-Minlo likewise reproduce NLO accuracy for Bj quantities, we are interested to resum the k t -jet resolution variables y 01 and y 12 , which we here define as where Q B and Q BJ are, in the context of the Caesar resummation framework [55], the hard scales of the problem, largely determined by the respective Born kinematics. We take Q 2 B = m 2 B , where m B is the invariant mass of the system B, and Q 2 BJ = Q 2 B v 01 = y 01 . With this definition the resummation of y 12 amounts to -up to corrections owing to the lack of monotonicity of the clustering sequence in d (n) , which we neglect -a resummation of large logarithms of the transverse momentum of the second hardest relative to the hardest branching in the event. For what follows we notate these large logarithms In discussions and formulae that apply equally well to the Bj-Minlo and Bjj-Minlo computations we simply use L to refer to either L 01 or L 12 . Equally, we will use y to ambiguously mean y 01 and y 12 , and v to mean v 01 and v 12 , when safe to do so. We now introduce the kinematic variables specifying the hard configurations about which we intend Minlo to resum the y 01 and y 12 variables. First consider applying the k t -jet algorithm to events such that they are clustered to the point of containing just a single jet (pseudoparton) and the system/particle B. We define directly from such ensembles Bj underlying Born variables, where the setΦ B specifies the configuration of B in its rest frame, including its invariant mass, 2 y J is the rapidity of the jet, y B is the rapidity of B, p J T the transverse momentum of the jet, and φ J its azimuthal angle. After a subsequent clustering with the k t -jet algorithm the jet/pseudoparton is also removed leaving just the system B, for which we further The definitions of Φ BJ and Φ B can also be considered as projections from real (or multiple emission) kinematics onto Born kinematics for Bj and B final-states respectively; note that, strictly speaking, in that context the jet in the projected Bj kinematics should be understood as being massless. The choices of y B and y J are motivated by our expectation that even basic formulations of the Bj-and Bjj-Minlo calculations will reproduce well the shapes of these quantities, as they are predicted in the respective (conventional) NLO B and Bj computations. Our choosing of p J T in Φ BJ is made in light of the fact that this variable, as defined here, is equal to √ y 01 , which we expect to greatly increase the consistency between the Bj-Minlo (which resums precisely v 01 ) and the Bjj-Minlo calculations, which we intend to have identically reproduce NLO Bj predictions such as p J T , as defined here, according to the exclusive k t -jet clustering algorithm. The latter consideration is important in the context of nesting the v 01 and v 12 resummations, with a view to having a Bjj-Minlo simulation NLO accurate for Bjj, Bj, and inclusive B production processes.

NNLL σ resummation
By differentiating and expanding the master resummation formula in ref. [55], we are able to derive simultaneously LL and NNLL σ accurate expressions for the v 01 and v 12 spectra in B and Bj production processes respectively. In a nutshell, one takes the NLL resummation of ref. [55], matched to NLO, and proceeds to omit NLL terms O(ᾱ 3 S ) and beyond in the resummed exponent, specifically those due to observable-dependent multiple emission effects. Details on these manipulations can be found in appendix A.1. The general expression we derive can be written simply in the form [56]: 3 where L is our luminosity factor .

(2.2)
Since this form applies to both jet resolutions in B and Bj production, the components inside it should be understood as referring to one of these two processes, e.g. Φ and L refer to Φ BJ and L 12 in the Bj case, and Φ B and L 01 in B-production. Equally, v refers to v 01 in the latter case and v 12 in the former.
First let's overview the resummation formula, eq. 2.1, before disappearing into the details. The first factor in eq. 2.1, dσ 0 /dΦ, denotes the leading order cross section for B or Bj processes as appropriate. Within dσ 0 /dΦ the scale used for the evaluation of the parton distribution functions is µ F and the scale in any implicit strong coupling constant factors is µ R . The function H 1 includes hard virtual corrections to dσ 0 /dΦ. The Sudakov form factor is present in eq. 2.1 as exp [−R (v)]; here we have made a simplification with respect to the notation of ref. [55], including in its definition contributions from soft-wide angle radiation and observable-dependent multiple emission corrections.
The q ( ) x , µ 2 terms in the luminosity factor, eq. 2.2, are parton distribution functions (PDFs), for a given incoming leg, , with momentum fraction x , evaluated at scale µ. The product of PDF ratios runs over n i = 2 incoming legs. The functions C 1 involved in convolutions with PDFs in the luminosity factor, L, are due to universal hard-collinear corrections. The renormalization and factorization scales µ R and µ F are understood as being ∼ Q [55].
To try to lighten the formulae we use the following abbreviations, The Sudakov form factor exponent in the NNLL σ differential cross section formula (eq. 2.1) is given by The G ij contributions are due to independent soft-collinear / collinear emission contributions, they are given by where for a quark leg , The sum, , runs over all n hard colour-charged legsn = 2 for B-production, n = 3 for Bj.
Single logarithmic soft-wide angle emission contributions are included via the S 1 term. Soft-wide angle radiation is obviously sensitive to the structure of the underlying hard event on large angular scales, so in contrast to the collinear contributions above, this piece is sensitive to the orientation of the hard external legs and not just their charges. For B (n = 2) and Bj(n = 3) processes we have where Q ij = |2p i .p j |. For the case of two/three hard gluon legs, we simply replace C F by C A in S 1 and, in addition, q, q , g with g 1 , g 2 , g 3 (see bottom of pg. 38 in ref. [55]). By writing S 1 in this form for the n = 3 case one can already glimpse, in the first term, its interpretation in terms of coherent emission from the n = 2 kinematic underlying the n = 3 one; we discuss this in more depth later on.
The K in the O ᾱ 2 S part of R is the two-loop cusp anomalous dimension Concerning the Sudakov form factor, the only remaining part needing introduction is F 2 . The F (R ) function of ref. [55] accounts for NLL corrections arising from a resummed observable's sensitivity to multiple emission effects; for observables whose behavior is largely dictated by the leading single emission F (R ) → 1. The Caesar F (R ) factor is understood to depend only on the flavours of the particles entering and exiting the hard scattering and the multiple emission properties of the observable, it does not depend on the kinematics of the underlying hard scattering (Φ) [57].
To NLL/NNLL σ accuracy F (R ) = 1 for the v 01 resummation (B-production). The combination of factors F 2 (2G 12 ) 2ᾱ 2 S L 2 is the next-to-leading term 4 in the fixed order expansion of the F (R ) function and as such defines F 2 . From refs. [38,57,58] we derive the following process-independent expression for F 2 , for jet rates in the exclusive k t algorithm (2.10) We have tested this expression using the numerical implementation of the Caesar formalism for resummation of y 23 in hadronic jet production and y 12 in hadronic Z boson production. With the exception of the qg and gq channels in Z production, for which only 3% differences were found, our F 2 expression yielded agreement with the Caesar program at the per mille level in all tested processes and channels.
In the resummation formula, eqs. 2.2-2.1, for the PDF dependent pieces we have adopted the notation q x, µ 2 where P ij (x) are the regularized leading order Altarelli-Parisi splitting functions (see e.g. appendix A.3 of ref. [55]). We also identify q ( ) x , µ 2 = q ( ) i x , µ 2 , with i the flavour of the hard parton with momentum p , and we employ the following notation to denote matrix multiplication and convolution in x-space The last things we need to introduce in our resummation formula, eq. 2.1-2.2, are the H 1 and C 1 terms. To this end we first define the cumulant, Σ R , of the NNLL σ resummed spectrum as Since dσ R is expressed as a total derivative we quickly find the following approximation to the NLO B/Bj production cross section: The cross section dΣ R,1 /dΦ| H whereχ 1 (Φ), 5 being localized at v = 0, encodes the regular virtual and the hard collinear contributions, with dσ F ,1 being the real emission contribution to dΣ NLO (L), with its v → 0 end-point subtracted and included inχ 1 . From eq. 2.15 we obtain directlȳ We separate hard-virtual and hard-collinear corrections inχ 1 as follows: where the C 1 terms represent the contribution due to hard-collinear splitting in the initial-state and H 1 is the remainder, including the hard-virtual component. H 1 contains terms canceling the µ R dependence, while C 1 has terms which correspondingly compensate the µ F dependence, of dσ 0 /dΦ B .
The precise details of these terms are irrelevant for the implementation of the method being proposed (this can be considered one of its advantages), so we can safely leave further specification of them to appendix A.1. We only stress that in the C 1 function, in eq. 2.2, which is convoluted with a PDF evaluated at scale µ F √ v, the explicit factorization scale is µ F , not µ F √ v, i.e. C 1 in eq. 2.2 is precisely as it is written in eq. A.13, so the derivative with respect to L in eq. 2.1 passes through C 1 , only acting on whatever follows it.
The structure of eq. 2.1 with regard to the inclusion of theχ 1 term, can be intuitively understood by considering that the Sudakov and PDF factors in eq. 2.1 are resumming the effects of all orders 5 What we have denotedχ 1 (Φ) ref. [59] denotes as C soft/collinear radiation around the hard scattering, described by dσ 0 in the case of the leading term, and, that the same pattern of radiation also occurs with respect to the hard scattering including hard radiative effects dσ 0χ1 . In other words, one can view the resummation as being taken with respect to essentially two separate processes and then adding these two resummations together, one process being the higher order analogue of the other. The fact that the resummation should be identical with respect to either process can be understood by considering that the soft longwavelength radiative corrections -encoded by the Sudakov form factor, running coupling and PDFs -will not be able to probe the internal details of the hard scatterings they attach to.

Comparison to Ckkw/Minlo
Before continuing, it is worth comparing the formulae presented here to those of Ckkw/Minlo [14,26], in particular the Sudakov form factors. The Sudakov form factor exponent in the latter articles 6 for a collection of (pseudo-)partons (indexed by ) evolving from some scale Q down to a resolution scale y, is given by The integral on the right-hand side of eq. 2.18 is the difference between the total Sudakov form factor exponents used in the Ckkw/Minlo prescription and that proposed here based on Caesar For B production F 2 = 0 and S 1 ∝ ln m B Q B , with Q B set equal to m B in the original Minlo proposal, hence, in this case, the second term on the right of eq. 2.18 vanishes. Thus, in B production the Sudakov form factor of the original Minlo procedure is fully consistent with that prescribed by Caesar.
For the Bj case, not forgetting that here in section 2 we are restricting ourselves to considering the region y 01 O(m 2 B ), F 2 is not zero, and S 1 has non-trivial dependence on the underlying Bj kinematics. Therefore, in the region where our Caesar-based formula is strictly valid we have a discrepancy between what is suggested by it and by Minlo. In particular the original Minlo proposal has omitted NNLL σ terms due to multiple emission corrections (F 2 ) and, more importantly, NLL σ contributions due to soft-wide-angle radiation (S 1 ). Thus, in the region y 01 O(m 2 B ) Bjj-Minlo, implemented according to the original proposal in ref. [26], would formally not be LO accurate in the description of Bj-inclusive quantities, with ambiguities arising between it and conventional LO of order √ᾱ S times the leading order term. With the benefit of hindsight it is perhaps obvious that the original Minlo procedure would have this problem in this region, since we know that its Sudakov form factors contain only soft-collinear and collinear terms, yet soft-wideangle radiation from a Bj state will be logarithmically enhanced too, even if the underlying Born partons are widely separated.
In section 3 we also consider this comparison (for Bjj-Minlo) in the region y 01 < O(m 2 B ).

Minlo jet resolution spectra
In the Minlo framework, in all cases, we start with an NLO cross section: for the v 01 resummation in B-production our fundamental ingredient is the NLO Bj cross section, while for v 12 resummation in Bj-production it is that for NLO Bjj. We write these cross sections as a sum of a part which is finite as v → 0, dσ F , plus a singular part obtained by expanding the NNLL σ resummation formula (eq. 2.1) dσ S , and a further singular-remainder piece, dσ SR , which is defined as all singular terms which were not already contained in dσ S : Expanding the resummed differential cross section up to and including O ᾱ 2 S terms, we obtain where the explicit H nm coefficients are documented in the appendix A.2. Since the resummation formula we used to derive this fixed order expansion was NNLL σ accurate, it only predicts part of the full N 3 LL σ coefficient, ∼ᾱ 2 S , thus we have a singular remainder term, where dσ R is the resummation cross section, eq. 2.1, a total derivative, and dσ MR holds all remaining large logs: (2.23) In eq. 2.23 the K R/F ∈ 1 2 , 2 denote rescaling factors applied to the renormalization and factorization scales, µ R/F , defined at the start of the Minlo procedure (see A.3 for details), for the purposes of assessing scale uncertainties. The last term in eq. 2.22, dσ F , is more precisely dσ MF , the replacement dσ MF → dσ F being made on the grounds that the Minlo operations preserve the fixed order expansion up to and including NLO terms, as well as the fact that dσ F (and dσ MF ) is finite for v → 0.
Since R 21 = 0, the Minlo jet resolution spectra in eqs. 2.22 are equal to the NNLL σ jet resolution spectrum in sect. 2.2 (dσ R ) up to N 3 LL σ differences.

Integrated Minlo jet resolution spectra
Making use of the fact that dσ R is a total derivative with respect to L (eq. 2.22), and the definitions ofχ in terms of H 1 and C 1 , it is fairly straightforward to show 7 that on integrating over all v The contaminating dσ MR term consists of a N 2 LL σ piece, ∝ R 21 , and N 3 LL piece ∝ R 20 −β 0 H 1 .
For the regions in which the Caesar formalism holds R 21 = 0, as discussed under eq. 2.21.
If we assume that we were ignorant of the value of R 21 , dropping terms over which we have no control, i.e. beyond NNLL σ order, we can neglect the L dependence ofᾱ S and PDFs in dσ MR , and all but the leading term in the Sudakov form factor exponent ∝ G 12 . With these approximations the dσ MR integral becomes: (2.25) The O(ᾱ 3/2 S ) ambiguity in eq. 2.25 attributes to neglect of N 3 LL σ terms. So, if our knowledge of R 21 be wrong, for whatever reason, the Minlo inclusive cross section would deviate from the exact NLO one by terms of order O (ᾱ S ) relative to the LO contribution (dσ 0 ).
Sticking to the regions for which the Caesar formalism holds, our starting resummation formula and the Minlo cross section formulated with it is NNLL σ accurate, i.e. R 21 = 0 and our ignorance is located downstream in the N 3 LL σ terms ∝ R 20 −β 0 H 1 . Dropping terms now only of N 4 LL σ accuracy we can again neglect the L dependence of the coupling constant and PDF terms, and all but the leading double log term in the Sudakov form factor, giving All large logarithms in this modified Minlo spectrum, eq. 2.28, are wrapped up as a total derivative and it's trivial to verify (more-or-less exactly as in appendix A.4) that the integral over all v (L) gives the conventional NLO cross section without any spurious terms, as were examined in sect. 2.5.
Thus we interpret the spurious terms that arise on integration in sect. 2.5, as being due to neglect of any NNLL σ (for the scenario R 21 = 0) and N 3 LL σ terms in our Minlo Sudakov form factor (eq. 2.4). To remove the spurious terms and recover NLO accuracy on integration over L we should try and include these terms in the latter. The Minlo approach of ref. [28] with −∆R (v) exactly as in eq. 2.27, leading to The modification to the Sudakov form factor in eq. 2.27 is equal to that in eq. 2.29 with differences only starting at the N 4 LL σ level. Thus, in this modified Minlo spectrum, eq. 2.30, the integral over all v (L) gives the conventional NLO cross section without any spurious terms.

Alternative approach to Minlo
The message from eqs. 2.27-2.28 and eqs. 2.29-2.30 is the same: including the appropriate corrective factor on top of the default Minlo Sudakov form factor, exp [−R (v)], we recover from Bnj-Minlo NLO accuracy also for Bmj inclusive observables (m = n − 1). We now suggest to turn around the latter fact and use unitarity to effectively determine the missing piece of the Sudakov factor, exp [−∆R (v)], at a level of accuracy sufficient for our aims.
Minded by the equivalence between the Minlo formulations in eqs. 2.27-2.28 and eqs. 2.29-2.30, we describe now how to implement, approximately, the 1 − ∆R (v) factor of the latter, without explicit knowledge of the R 21 , R 20 and H 1 terms. We define what we consider to be the discrepancy in the Sudakov form factor at NNLL σ as: In eqs. 2.31-2.32 we abbreviatedᾱ S Q 2 →ᾱ S . The leading term in the integrand of the denomi- The choice of the h (L) function in eq. 2.31 is not a rigid one, as we comment on later. The freezing parameter ρ in eq. 2.32 is taken to be ∼1 by default. From eqs. 2.24-2.25 we get an expression for the numerator of δ(Φ), and by similar approximations to those used for the latter, an expression for the denominator (appendix A.5), giving overall Since δ(Φ) is an order one quantity (provided ρ 1) we can safely define the following modification of the original Minlo distribution i.e. the corrected Minlo distribution precisely returns the true NLO inclusive cross section on integrating out the radiation, unambiguously. This correction is achieved while leaving the NLO accuracy of the input cross section intact; the weighting factor in square brackets in eq. 2.34 being The modification in eq. 2.34 also does not interfere with the Minlo cross section at NLL σ . 8 The ρ parameter guards against the 1 − ∆R (v) approx factor in eq. 2.34 becoming negative, which can happen in the regionᾱ S L 1, if δ (Φ) is positive, leading to an unphysical spectrum at v → 0.
We remind that the regionᾱ S L 1 has been anyway, from the beginning, outside the control of our calculational setup, which is only adequate down to the regionᾱ S L 2 ∼ 1. Introducing ρ also tames the integrand in the denominator of δ (Φ), so it can be determined/applied by simply weighting events appropriately in analysis of the original dσ M distribution, without issues of numerical convergence.
To provide some advance reassurance, should any be needed, in our feasibility study in section 4 we carry out what we consider to be broad variations of the ρ parameter, finding our results 8 To see this consider re-expressing the square bracket term in eq. 2.34 as an exponential.
exhibit marginal sensitivity to it, in regions of practical interest.
In the region of applicability of Caesar, our knowledge of the spectrum is complete at NNLL σ , i.e. we know that if the Sudakov form factor in eq. 2.4 is implemented in Minlo, R 21 = 0. All equations and analysis above remain valid for R 21 = 0 though. The latter implies δ (Φ) (eq. 2.33) (1), meaning the correction factor, eq. 2.34, has effectively less work to do. For R 21 = 0 the latter correction factor clearly affects the spectrum at NNLL σ . However, we argue that since, for R 21 = 0, the coefficient of the correction in eq. 2.
for the domain of validity of the Caesar formalism, it may have seemed more natural to have made h (L) ∼ᾱ 2 S L in our discussion here, instead of ∼ᾱ 2 S L 2 , however, for the formal reasons just discussed, we see no great advantage in doing so. Furthermore, we plan to employ the method also in the region where Caesar is not valid, in the next section, so having a more widely applicable h (L) function, which nominally assumes the distribution it is correcting is NLL σ accurate, is preferable.
It may be tempting to think that one can also apply this procedure even if the initial input Minlo distribution was only LL σ accurate, supplying, in that case, the h (L) function with one more power of L, in order to keep δ (Φ) ∼ O (1). While this appears compatible with the recovery of NLO Born kinematics, maintaining also NLO accuracy of the initial Minlo simulation, the expansion of the product of the latter factor and the initial LL σ Sudakov form factor has a different functional form to that of a NLL σ Sudakov form factor. In other words, one cannot then view the resulting correction (eq. 2.34) as approximating missing higher order pieces of the Sudakov form factor, which was has been our guiding principle throughout. This conflict can only be resolved by making h (L) ∼ᾱ S L, however, in that case the correction factor (eq. 2.34) will clearly violate NLO accuracy of the initial Minlo program. We therefore consider it a requirement that the relevant resummation in the initial, uncorrected, Minlo program be at least NLL σ . Fortunately, this is a rather low theoretical threshold to cross by today's standards, and really the only non-trivial NLL σ ingredients required are the (soft-wide-angle) S 1 Sudakov coefficients (eqs. 2.4-2.7). It is well understood how to obtain the latter soft-wide-angle pieces, and it is not a particularly onerous task to do so nowadays. Indeed there is much publicly available, automated, machinery which can be straightforwardly adapted to this end, e.g. in Powheg-Box [4] and Madgraph5_aMC@NLO [7].
Finally, it is also the case that the aforementioned S 1 terms are trivial for processes where the underlying Born comprises only two or three coloured partons (e.g. Bj-and Bjj-Minlo).
There are numerous possible variations, tangents and refinements one can explore along the lines presented here, all leading to eq. 2.36, with or without ambiguities. For example, one can easily enough conceive of modifications which avoid the introduction of the parameter ρ. Equally, there are other ways to view the formulae in this section, most of which are obvious. We do not want to digress, to avoid diluting the basic idea and straying too far from the goals in the introduction.
In particular, we choose not to discuss to what extent we have formally improved the description of the resummation region, but rather we now get on with demonstrating the practicality of the above, and its extension beyond the merging of two units of multiplicity.

Merging three units of multiplicity
We now turn to address the problem of getting the Bjj-Minlo calculation to return NLO predictions for inclusive B-production observables, as well as Bj and Bjj inclusive quantities.
Since Bj-Minlo contains at most one final state parton with NLO accuracy, from now on we discuss modifications to Bjj-Minlo only. Furthermore, we focus on modifications needed to address the remaining problematic region y 12 y 01 m 2 B , with y 01 m 2 B having been covered in sect. 2.
Where necessary, we use a superscript [01]/ [12] on quantities G ij , S 1 etc to distinguish those associated to the y 01 resummation, from those associated to that of y 12 .
The basic idea here is simple and can easily be improved; we briefly discuss such refinements later on, with a view to future work. We require that a lower multiplicity Bj-Minlo /Nnlops simulation has already been built, e.g. with the procedure of sect. 2.6, or along the lines of ref. [28], recovering NLO accuracy for Bj-and (NNLO) B-inclusive observables. We propose to apply method of sect. 2.6 to the Bjj-Minlo simulation, with the obvious replacement dσ M → dσ BJJ M therein, but also with the conventional fixed order distribution dσ NLO replaced by dσ BJ M . It then follows, trivially, that Thus, the resulting Bjj-Minlo distribution is targeted onto the Bj-Minlo inclusive Φ BJ distribution, without diminishing its own NLO accuracy. In this way the Bjj-Minlo simulation can be made NLO accurate for Bj-and (NNLO) B-inclusive observables.
The essential point one needs to prove for the self-consistency of the method in this context is the same one as in sect. 2.6, i.e. that the δ(Φ BJ ) that gets extracted does not blow up and risk breaking the NLO accuracy of the initial uncorrected Bjj-Minlo. This basically boils down to saying that the existing Minlo procedure resums v 12 and v 01 both with NLL σ accuracy. As discussed at the end of sect. 2.6, if all we cared about was unitarizing the cross section, this requirement could be loosened to that of having just LL σ accuracy in place, including a further power of L in h (L). The price of that ignorance would be that the correction can no longer be interpreted as an approximation to missing higher order contributions in the Sudakov form factor, i.e. one essentially gives up on a physical interpretation of the mechanism of unitarity violation and, correspondingly, one begins to warp the spectrum by higher order ambiguities that bear no relation to any kind of resummation.
However, as we go on to explain, we understand the Bjj-Minlo cross section meets already the above NLL σ specification, with the exception of a sub-leading kinematic region, which should not be difficult to accommodate.

NLL σ resummation
It is a general underlying assumption of the Caesar formalism that the Born configurations (Φ BJ in this case) consist of hard, well-separated, partons. So the NLL/NNLL σ theoretical framework from which we derived the resummation formula, that was the starting point for section 2, is not guaranteed to hold here, where we also need control v 12 resummation in the region √ y 01 m B .
In section 3.1.1 we argue that the Caesar resummation formula for v 12 = y 12 /y 01 resums large logarithmsᾱ n S L m 12 , m ≥ 2n − 1, independently of the value of v 01 = y 01 /m 2 B , i.e. even in the region y 01 m 2 B . This is based on the following two considerations: i.) it is straightforward to show that the Caesar Sudakov form factor for v 12 resummation is equivalent to that prescribed by the coherent parton branching formalism at NLL σ , except for a sub-dominant subset of soft wide-angle radiation contributions, beyond the accuracy of the latter formalism; 9 ii.) the leading NLL σ terms in the expansion of the Caesar v 12 cumulant distribution, are determined by integrating the radiation pattern of a single soft/collinear emission relative to an emitting BJ state, over the ordered region y 12 < y 01 , and this pattern/integral derives independently of whether y 01 m 2 B or not.
In i.) we are saying that the Caesar Sudakov form factor must be at least LL accurate for v 12 , since it agrees with the analogous expression derived from the coherent branching formalism, which is understood to have at least that accuracy, regardless of the value of the underlying Φ BJ configuration. 10 Accepting i.) and ii.) together then implies that the Caesar resummation formula is NLL σ regardless of the value of the underlying Φ BJ configuration: since the Sudakov form factor is present in the resummation formula as an overall factor, if the expansion of the formula generates just the leading NLL σ terms in the cross section correctly, it generates all of them correctly. For the reader who is willing to accept the statements above without detailed explanation (the first of which is not obvious) we recommend skipping 3.1.1.
In section 3.1.2 we go on to include v 01 resummation at NLL σ . To this end we notice how, if we include on top of the Caesar v 12 resummation formula, matched to leading order Bjj, also the v 01 Sudakov form factor, on integrating out y 12 we obtain the Bj-Minlo distribution to NLL σ .
The analysis in sect.

has made it clear already that this Bj-Minlo distribution recovers the
Caesar v 01 resummation formula on further integration over the rapidity, y J , and azimuth, φ J , of the remaining pseudoparton.
With the latter modification we come full-circle: by concatenating the two Caesar resummations we get the same resummation as the original Ckkw [14,27] and Minlo articles [26], to NLL σ , modulo the terms in the v 12 Sudakov form factor mentioned overhead in item i. If we restrict ourselves to the same accuracy remit as the coherent branching formalism aims at, the only part of our prescription not already specified in the original Ckkw paper [14], is the inclusion of the PDFs.
Again, our argument to extend the prescription to include PDFs (and also the aforementioned wide angle terms) is based on the idea that if the resummation formula carries an overall LL accurate Sudakov form factor and reproduces just the leading NLL σ terms correctly, it surely reproduces all of the NLL σ terms in the cross section. It is reassuring then that our extension, in that respect, to take into account the PDF effects, is also consistent with that of the Ckkw paper on hadronic collisions [27], and the original Minlo prescription [26].
While the arguments behind our nested resummation are strong enough to convince us of its correctness, we do not consider that we have definitively proven it.
Since it can lead to confusion, in reading this subsection 3.1, the sole goal of which is to give a conjectured NLL σ resummation formula, we advise the reader to temporarily abandon all thoughts about matching to NLO, and imagine instead just matching to LO Bjj cross sections.
For what concerns this subsection we focus on large L 12 logarithms relative to a given Bj state. Large L 01 logarithms are discussed next in sect. 3.1.2. The statements we make regarding L 12 resummation should hold independently of the behaviour of the dσ/dΦ BJ underlying Born cross section, be it pathological or otherwise. However, should reassurance be needed already, the divergent y 01 → 0 behaviour of dσ/dΦ BJ is ultimately tamed by inclusion of a Sudakov form factor consistent with Caesar and the coherent parton branching formalism. We now specify the connection between the Sudakov form factors in the Caesar approach and those used for the coherent parton branching formalism/Ckkw [60][61][62][63][64][65][66][67][68][69][70][71]. The latter formalism is capable of resumming v 12 logarithms also in the small v 01 region. The key point that comes out of this analysis, in regards to making the case for the nested resummation, is that (ignoring potentially enhanced ln z terms 11 ) the Sudakov form factors associated with the y 12 resummation in both approaches are the same to NLL σ .
For processes with n = 3 hard legs, all S [12] 1 coefficients (eq. 2.8) can be written, without approximations, as a piece containing a logarithm of y 01 plus a remainder term, ∆S 1 , which, crucially, for Q B = m B , Q BJ = √ y 01 , has no large-logarithmic dependence on y 01 : Rewriting S [12] 1 as in eq. 3.2 is the key to understanding the connection between Ckkw and Caesar here. In eq. 3.2 the G 12 coefficient is that which one would write down for the n = 2 process underlying the n = 3 one; qq → W/Z and gg → H for jet-associated W/Z and Higgs boson production processes. Explicit expressions for ∆S 1 are given in appendix A.6.
While ∆S 1 is free of large y 01 logarithms, it is not zero. In the n = 3, 2 → 2, hard configurations with a gluon in the final-state, ∆S 1 contains terms proportional to ln z, where z = m 2 B /ŝ, withŝ the invariant mass of the 2 → 2 collision. In the n = 3, 2 → 2, hard configurations with a fermion emitted in the final-state, also terms proportional to ln (1 − z) are present in ∆S 1 . Such terms are thrown out in the coherent parton branching formalism as being beyond the accuracy aimed at there for exclusive quantities, namely, control of all termsᾱ n S L p 01 L q 12 , p + q ≥ 2n − 1; heuristically, that accuracy implies a resummation of an infinite number of soft and collinear emissions with, in addition, up to one soft-wide-angle, or hard-collinear emission. Thus, in order for soft-wide-angle emissions, which the S [12] 1 terms are to account for, to be within the accuracy remit they must have been emitted from an underlying Bj state for which z → 1. Equally, in the case of the n = 3 reactions with a fermion in final-state of the underlying Born, the fact that the n = 3 state is arrived at by fermion emission, means that further radiation must be soft-collinear to register within the formalism's accuracy. Thus the ∆S 1 terms are (so far) outside the scope of the coherent parton branching framework, indeed, they would appear to exactly the kind of "large angle soft gluon 11 The Caesar framework (like many other works) neglects the potential small x problems anyway.
contributions of order α n S L 2n−2 " which the formalism neglects in the case of multi-jet distributions (pg. 11, ref. [65]).
In comparing the Caesar formulas to those of the coherent parton branching framework we therefore now drop ∆S 1 terms, and those contributing beyond NLL σ , in the Sudakov form factors, leading to the replacements Translating eqs. 3.3-3.4 in terms of the notation of the coherent parton branching formalism we get, without approximations, (3.5) where ∈ [01] means one of the two coloured legs which directly attaches itself to B, 12 while ∈ [12] means any of the three coloured legs external to the Bj state. Definitions of the Sudakov form factors ∆ are given in appendix A.7, they are the same as those used widely in the literature on the coherent parton branching formalism/Ckkw (e.g. ref. [14]). Observe how the form of the product of the two Caesar-style Sudakov form factors gives the breakdown one expects in terms of the Sudakov form factors employed by the coherent parton branching formalism/Ckkw method.
Continuing to neglect ∆S 1 terms, the y 12 Sudakov form factor can be rewritten without further approximation as Finally, we have that the coherent parton branching formalism and the Caesar y 12 resummation formula are consistent in regards to the Sudakov form factors they would assign for the y 12 (and y 01 ) resummations, at the level to which the former is accurate. Caesar's accounting for soft-wide angle resummation, via the leading part of its S [12] 1 term (eq. 3.2) is essential for this non-trivial agreement. Beyond the domain of validity of the coherent branching formalism we only have LL agreement with the coherent parton branching formalism in the Caesar Sudakov form factor for y 12 resummation.
Following the argument laid out surrounding bullets i.) and ii.) in sect. 3.1, we therefore consider the Caesar y 12 resummation formula to be NLL σ accurate also in the region y 01 m 2 B . For this to be false requires either: i.) the statement that the y 12 resummation formula is not LL accurate for arbitrary Φ BJ to be false, which conflicts with the coherent parton branching formalism; ii.) the leading NLL σ terms in the expansion of the Caesar y 12 cumulant in the region y 01 m 2 B do not follow directly from integrating the soft/collinear radiation pattern of a single emission with respect to the emitting Bj configuration, over the region y 12 < y 01 .

y 01 resummation
Going back to our initial resummation formula of section 2, neglecting higher order terms, we now understand the following resummation formula to be NLL σ independently of Φ BJ , [12] , y 12 q ( ) x [12] , y 01 , where x [12] refers to the momentum fractions of the incoming partons colliding to make the Bj system. Integrating this formula over y 12 we obtain the leading order Φ BJ distribution, dσ BJ 0 /dΦ BJ , up to NLO-sized ambiguities. The renormalization scale in the coupling constants in dσ BJ 0 is µ R and the factorization scale in the PDFs is µ F . It then follows directly (given the correspondence between the Minlo procedure in sect. 2.4 and the initial resummation formula eq. 2.1) that for we reproduce the Bj-Minlo distribution, and hence also the Caesar y 01 resummation formula, to NLL σ accuracy. We conclude that dσ BJJ R in eq. 3.8 above, is NLL σ accurate in the resummation of L 12 for arbitrary given Φ BJ , and that it reproduces, on integration, the L 01 resummation to the same precision.

Bjj-Minlo jet resolution spectra
Expanding the conjectured resummation formula inᾱ S to give the associated NLO approximation for the Bjj cross section, we get where the coefficients H [12] nm have the same form as those introduced in sect. 2.4, with the renormalization and factorization scales µ R = m B and µ F = √ y 01 throughout; 13 explicit expressions for these can be found in appendix A.8. If we now trace the effects of the Minlo procedure on this fixed order expansion, by analogy to the exercise of sect. 2.4, we find that the final Minlo cross section is in agreement with the resummation formula, eq. 3.8, up to sub-leading terms outside the control of the latter. Specifically, here the Minlo procedure for the NLO Bjj cross section is: 13 Including inside the PDF factors of the dσ BJ 0 term.
3. Multiply by the Minlo Sudakov form factors andᾱ S ratios: With these operations we find we can write dσ BJJ M = dσ BJJ R , as in eq. 3.8, neglecting sub-leading terms unaccounted for by the dσ BJJ R formula.
Recalling that the product of the Caesar Sudakov form factors is equivalent at NLL σ accuracy to the product of those prescribed in the Ckkw method and the original Minlo procedure [26], modulo the ∆S Sudakov form factor terms. As indicated already in the introduction to this section, the 'product' of the two Caesar resummations has returned us, somewhat remarkably, almost exactly back to the Ckkw/Minlo recipe.

Integrated Bjj-Minlo jet resolution spectra
Granted that the Bjj-Minlo procedure is NLL σ accurate for the v 12 resummation and NLL σ for v 01 when y 12 is integrated out, it follows that where the e n coefficients are O (1). The e n coefficients carry no divergent 1/y 01 factors -these cancel between the numerator and denominator of δ 14 -equally, they contain no large L 01 factors. 15 Thus δ(Φ BJ ) is formally also O (1), neglecting deep Sudakov regions whereᾱ S L 2 01 1. This means that we are justified in applying the procedure of sect. 2.6 with the replacements dσ M → dσ BJJ M , dσ NLO → dσ BJ M , leading to eq. 3.1, and hence recover also NLO accuracy for B-inclusive quantities, or NNLO accuracy, should the Bj-Minlo distribution have been reweighted to NNLO [29].

Feasibility study
In the following we show how the above merging of three units of multiplicity works in a practical implementation. For this we consider Higgs production at the LHC, with a collision energy of 8 TeV. We 'merge' the Hjj-Minlo simulation to an existing Hj-Minlo simulation reweighted to NNLO according to the prescription of ref. [29]. In the following we therefore make predictions that are NNLO accurate for inclusive Higgs boson production, and NLO accurate for Hj and Hjj observables. The inclusive matrix element predictions are matched to the parton shower using the Powheg method.

Implementation
In order to simplify the implementation and require no changes to the existing Hjj and Hnnlops processes in the Powheg-Box, we have chosen to work at the level of the Les Houches events (LHE). The distributions formed from LHE in Powheg are NLO accurate, i.e. differences with a fixed order NLO computation are beyond NLO accuracy. Relatedly, they respect the NLL σ accuracy of Minlo: the difference in phase space w.r.t. the matrix elements is beyond LL and the matching to NLO is then enough to preserve distributions at the NLL σ level. We consider that working at the level of the LHE also simplifies the generation of the final results: we have written an independent code that reads in Hjj and Hnnlops LHE files and writes out a reweighted Hjj LHE file to achieve the results of the three units of multiplicity merging.
As described in sect. 2.6, we need to correct the dσ HJJ M dΦHJdL12 in such a way that when integrated over L 12 it returns the Hj-Minlo /Hnnlops Φ HJ distribution. This can be done by multiplying the fully differential Hjj calculation by the (1 − ∆R(v 12 ) approx ) factor as described in eq. 2.34.
This factor can only be computed after integration over L 12 , as is clear from eq. 2.31. To avoid performing the complete L 12 integration for every Φ HJ phase-space point, and given that this integral is too complicated to perform analytically, we instead have chosen to setup three three-dimensional interpolation grids for the three contributions to δ(Φ HJ ): the two terms in the numerator and the term in the denominator, respectively. The three dimensions are the rapidity of the Higgs boson, the rapidity of the hardest jet and the transverse momentum of the hardest jet. These being the dimensions making up the non-trivial part of the Hj phase space Φ HJ ; the dependence on the azimuthal angle, φ J , is completely flat. Indeed these interpolation grids can be filled quickly with the LHE from the existing Hjj and Hnnlops implementations in the Powheg-Box framework. We have generated the interpolation grids using rigid binning as well as a method based upon Parni [72] to dynamically create hypercubes in the three dimensions; we did not see appreciable improvements using the more involved Parni method and the results we present here are therefore based on the implementation using the simpler fixed interpolation grid bins.
In implementing the h (L) function of eq. 2.32, we have softened the abrupt transition at the freezing scale, manifested by the step functions. Specifically we implement h (L) as taking γ = 5. As γ → ∞ the h (L) function becomes exactly that of eq. 2.32, but for a rescaling ρ → ρ/(2 |G 12 |). Thus, h (L) becomes frozen in the region where the leading double log term in the Sudakov exponent is ≈ ρ. We probe the sensitivity of our results to ρ (and therefore, indirectly, also γ) by computing predictions with ρ = 1, , is beyond the scope of the coherent parton branching formalism, because the first emission, i.e. the one entering y 01 , is not enhanced by any large logarithm in that case. As clarified in sect. 2.3, it follows that Ckkw/Minlo does not lead to the correct Sudakov factors at NLL σ accuracy in the y 12 variable, in this region of phasespace: they miss the S 1 contribution due to soft-wide-angle radiation. In this feasibility study we chose to ignore these facts, with the expectation that technical issues might well instead present us with more serious, immovable obstacles. Formally, if these missing contributions would turn out to be important, δ(Φ HJ ), with the definition of h(L 12 ) as in eq. 2.31, would no longer be an order one quantity, c.f. eq. 2.33. It is not difficult to include these terms in the Minlo framework, indeed one of the results of this work has been to pinpoint these and other terms, which can improve the quality of the resummation. We leave the implementation of such terms to future work, although all indications from the following results suggest that this stands to be, fortunately and unfortunately, a null exercise. We have checked that in all of phase-space δ(Φ HJ ) remains within the range of values associated to resummation constants used in Higgs boson transverse momentum resummation.
Lastly we comment that we work in the SM theory in which the top quark is integrated out.
This results in a well-known Higgs effective theory with tree-level interactions between the Higgs boson and gluons. This approximation breaks down if the Higgs or the gluons carry enough energy to resolve the shrunk top quark loop, e.g., when the leading jet transverse momentum exceeds the top quark mass. We also do not include b quark mass effects, setting the b quark mass and Yukawa coupling to zero. Accounting for these finite-quark mass effects at the Born level in the Hjj Powheg generator could be done in an analogous way to ref. [73].

Results and testing
In the hard matrix elements we set the Higgs mass to m H = 125 GeV and keep it stable throughout the simulation. The LHE are showered with Pythia6 [74], using the Perugia 0 tune [75] but with hadronization and multiple-parton interactions turned off. The central renormalisation and factorization scales are set according to the Minlo procedure. To assess the scale dependence in Hjj-Minlo we vary the renormalisation and factorization scales independently, by a factor two around their central values, omitting the two values where the scales are changed oppositely. This results in 7 curves, whose envelope gives the uncertainty band. For the Nnlops we use procedure advocated in ref. [29], resulting into 21 curves. In the merging of the new improved Hjj-Minlo results we keep the scales used in the input Hjj and Hj (which on its own is an input to Nnlops) calculations correlated. Hence, this also results in 21 curves, the envelope of which defines the uncertainty band. We employ the MSTW2008nnlo PDF set [76] for all contributions and refrain from showing uncertainties of PDF origin. We remind that the correction procedure, as described in sects. 2.6 and 3, should function such that quantities which are fully inclusive with respect to the y 12 variable have no sensitivity to ρ at all.
Thus, for at least the first few figures we look at in this section, focusing on fully inclusive and Hjinclusive observables, the aforementioned dark-red band should be (and is) invisible, being obscured by the horizontal black reference line. Moving on to more interesting observables, particularly probing the behaviour of the second jet/second pseudoparton in the event, the dark-red ρ-parameter band starts to emerge, but it is generally quite elusive.
We do not claim that variation of ρ, together with the renormalization and factorization scales, gives a realistic estimate of theoretical uncertainties in regions where large Sudakov logarithms occur. We content ourselves to say that ρ is an unphysical technical parameter introduced in our procedure, with systematics associated to it. We believe our variation of ρ, as described above, is a conservative estimate of these systematics, and we find them to be very much negligible.
Finally, statistical uncertainties are shown as vertical lines, however, for the most part these are negligible to the point of being invisible.

Inclusive quantities
In figure 1 we plot the rapidity of the Higgs boson; no cuts have been applied to the final state. The Hjj and Nnlops central predictions agree with one another to within 2%, with their uncertainty bands exhibiting a similar level of agreement. This indicates that the method and its implementation are performing as expected (eqs. 2.36-3.1). The uncorrected Hjj-Minlo prediction in blue is 10% away from the central Nnlops results, but this is fortuitous given that the scale uncertainty on the former is ∼ 30%. Moreover, given our theoretical analysis in the preceding sections of this paper, neglecting the sub-leading NLL σ ∆S 1 terms, we expect the Hjj-Minlo prediction here is only LO accurate, so the ∼ 30% uncertainty assigned to it is arguably too small. The uncertainty band associated to varying the ρ parameter as described at the beginning of this subsection 4.2 is so small that it is concealed within thickness of the black reference line in the upper right plot; indeed since this quantity is fully inclusive in L 12 , by construction of the procedure (sect. 2.6), the only way any such uncertainty could manifest here is as a result of technical problems and/or some statistical issues. In figure 2 we plot the Higgs boson transverse momentum spectrum. As with the Higgs boson rapidity distribution no cuts have been applied to the final state. Exceptionally, in this figure we compare Hjj and Hjj to the NNLL+NNLO predictions of the Hqt program [77][78][79][80][81], instead of Nnlops. Comparing Nnlops (not shown) and Hjj we find the two generators agree with one another to within 3% throughout the spectrum, except for the region p T 5 GeV, where the difference rises up to 15% in the p T < 2 GeV region. The latter differences owe to the finite size of the bins in our interpolation grids, coupled with the fact that the distribution is changing very rapidly for p T 5 GeV. Given this technicality, and the fact that this region is under poor theoretical control anyway, the conclusion, again, is that the method and its implementation work well. Turning then to the comparison with Hqt in figure 2, we see, pleasingly, that the method substantially corrects the shape of the pre-existing Hjj-Minlo simulation, with the resulting Hjj prediction agreeing very well with Hqt in the region where the latter is undeniably the superior calculation (p T 100 GeV). 16 In the high transverse momentum tail both Hjj and Hqt computations have the same NLO accuracy for this distribution. Differences between Hjj and Hqt occur there due to the different choice of scales in each code, roughly, p H T in the case of Hjj , compared to 1 2 m H in Hqt. The same comments made above for the Higgs boson rapidity distribution in regards to the uncertainty associated with the ρ parameter apply equally well again here.

Jet cross sections
In figure 3 we compare predictions for inclusive jet cross sections, between the Hjj (blue), Nnlops (dark green) and Hjj (red) generators, defined according to the anti-k t -jet algorithm [82] with radius parameter R = 0.4, for jet transverse momentum thresholds of 25, 50 and 100 GeV. In figure   4 we show the analogous set of plots for the corresponding exclusive jet cross sections. No rapidity cuts have been applied to the jets in making these plots. The behaviour seen in all of the inclusive jet cross sections of fig. 3, is as we would naively expect it to be. By construction, our improved Minlo method should reproduce Nnlops results essentially identically for 0-and 1-jet inclusive quantities (eqs. 2.36-3.1), while observables that receive their leading contributions from higher jet multiplicities are to be described as in the original Hjj-Minlo generator, which yields the more accurate predictions for those observables. Given that our Minlo improvement method is intended to return the 0-jet and 1-jet inclusive results of its 'target' Nnlops simulation, essentially without ambiguities, one might be tempted to ask why we can see even 2% differences between the Nnlops and Hjj predictions for the 1-jet inclusive cross sections. What the improvement procedure precisely does, without ambiguities, assuming a perfect implementation, is to have the improved Hjj result reproduce the Nnlops underlying Born kinematics Φ BJ (eqs. 2.36-3.1) which are defined by clustering events with the exclusive k t -jet algorithm, with R = 1.0. What is plotted in the 1-jet bins of fig. 3 is therefore not in one-to-one correspondence with the kinematics Φ BJ (consisting of a Higgs boson and a single pseudoparton in the final-state) but rather it is something which is also sensitive to additional radiation.
The Hjj and Nnlops generators are further in agreement as to the relative distribution of this additional radiation at the level ofᾱ 4 S terms, i.e. at the level of NLO corrections to Hj, however, at O(ᾱ 5 S ) differences do enter. Hence, even if the implementation were a perfect representation of our method, with infinite resolution in the Φ BJ grids, we can still expect to see differences between the Nnlops results and Hjj , for the 1-jet inclusive cross section, which are formally NNLO-sized in the context of the inclusive 1-jet calculation. This being the case, one can be quite satisfied with only 2% differences between the Nnlops and Hjj predictions for the 1-jet inclusive cross sections.
In fact, we examined the 0-and 1-jet inclusive cross sections, with a 25 GeV jet p T threshold, prior to interfacing with the parton shower, whereupon we found the 0-jet and 1-jet Hjj cross sections to be indistinguishable from their Nnlops counterparts, while the 2-and 3-jet bins remained identical to those of Hjj-Minlo. For the 0-jet exclusive cross sections in fig. 4 we see nice agreement between the Hjj and Nnlops predictions at the 1-2% level or better, as is to be expected by construction of our method.
To explain the 1-2% differences that can be seen we tender again the same theoretical explanation as above (the Nnlops and Hjj results differ, by construction, at the level of O(ᾱ 5 S ) terms), however, with such small differences we also cannot rule out imperfections in the implementation, e.g. artefacts due to the finite granularity of the grids and grid interpolation. We suffice to say that the differences between the Hjj and Nnlops computations of the 0-jet exclusive cross sections are negligibly small, while the unimproved Hjj-Minlo result sits 10-15% below them.
Lastly, we look to the the 1-jet exclusive cross sections. The plots in this case read that the Hjj prediction is different from the Nnlops one by 7% for the 25 GeV jet p T threshold, 5% for the 50 GeV threshold, and ∼ 0% for the 100 GeV threshold. Meanwhile, the unimproved Hjj-Minlo prediction is in agreement with the Nnlops prediction at the level of ∼ 0%, 10%, and 15%, for the same p T thresholds, respectively.
Since the Minlo improvement method we propose works to correct the inclusive 0-and 1-  σ (≥ 1 − jet) − σ (≥ 2 − jets). Clearly if σ (≥ 1 − jet) σ (≥ 2 − jets) differences in the latter will have limited impact on the exclusive 1-jet cross section. The latter scenario is enhanced by increasing the jet p T threshold and, sure enough, the pattern of the exclusive jet-cross sections seen in the case of the 100 GeV p T threshold, mirrors well what we see in the analogous inclusive jet cross section case, discussed overhead. To explain then the differences seen between the Nnlops and Hjj generators at the 25 GeV and 50 GeV jet p T thresholds, we note that the Nnlops exclusive 1-jet cross section is given by its inclusive 1-jet cross section minus its inclusive 2-jet cross section, on the other hand, by design, as can be verified in fig. 3, the Hjj exclusive 1-jet cross section is basically given by the Nnlops inclusive 1-jet cross section minus the Hjj-Minlo inclusive 2-jet cross section. Since the Hjj-Minlo inclusive 2-jet cross section at the 25 GeV jet p T thresholds is 10% lower than the Nnlops one, while the ratio of the inclusive 1-jet to 2-jet cross sections is roughly two, it follows that one can expect the Hjj 1-jet exclusive cross section to be 5% higher than the corresponding Nnlops one. Adding in the fact that the Hjj 1-jet inclusive cross section was already 1-2% above the corresponding Nnlops one, the 7% excess is actually very much in line with expectations based on how the method is intended to work, in particular, its preserving of the inclusive cross sections. A similar explanation holds for the 50 GeV jet p T threshold result, however, there the fact that the 2-jet inclusive cross section of Hjj-Minlo is low does not imply as big an increase is needed in the 1-jet exclusive bin to recover the 1-jet inclusive Nnlops result, since for that higher p T threshold σ (≥ 1 − jet) σ (≥ 2 − jets). We also remark that the Hjj exclusive 1-jet cross section results all agree with those of the Nnlops generator to within the thickness of the scale uncertainty bands.
We finally note that the dark-red scale uncertainty band associated to variation of the ρ parameter, as described in sects. 4.1-4.2, is invisible here: the effect of varying this parameter on these distributions is totally negligible.
To conclude the discussion on jet cross sections, we can say that all of the results we find are very much in line with expectations regarding how our method should function, all vindicating the method and its implementation.

Leading jet
In figures 5-9 we plot various quantities relating to the kinematics of the leading jet. In all of these figures jets have been defined according to the anti-k t clustering algorithm, with jet radius R = 0.4; no rapidity cuts have been applied to the jets. For figs. 7, 8, 9 jets were further defined as having a transverse momentum of at least 25 GeV. The results for the leading jet transverse momentum spectrum in fig. 5 read similarly to those reported for the Higgs boson transverse momentum spectrum ( fig. 2). The Nnlops and Hjj predictions agree very well throughout the spectrum, with the procedure correcting well for substantial (±15%) shape differences between the unimproved Hjj-Minlo result and the more accurate Nnlops prediction. Regarding differences between the Nnlops and Hjj results in the p T 5 GeV region, the explanation here is the same as for the case of the Higgs boson p T spectrum, namely, that the granularity in our discretized implementation of the Φ BJ phase space is not sufficiently fine to cope with the rapidly changing distribution for p T 5 GeV. We reiterate that this region is under limited theoretical control anyway. Indeed, rather than seek improved agreement of Nnlops and Hjj in the latter murky region, we might prefer to lessen the 3-5% deviation in the neighbourhood 60 ≤ p T ≤ 80 GeV. This region, where the Hjj-Minlo and Nnlops lines intersect, appears to be where the p T derivative of the difference between the two predictions is changing most rapidly, i.e. the numerator of δ (Φ BJ ) in eq. 2.31/3.11. It should therefore be possible to improve agreement between the Nnlops and Hjj results in this region by, for example, making use of (irregular) optimized grids and interpolation methods which can work on them. Overall, notwithstanding our unsophisticated implementation, agreement between the Nnlops and Hjj predictions is very satisfactory, providing significant improvement across the whole p T spectrum relative to the original  This is in contrast to the band associated with it in ref. [41], where additionally resummation scale and matching scheme variations were included in the envelope. Thus the JetVHeto error band here is considerably smaller than that shown in ref. [41]. We restricted the JetVHeto uncertainty estimate to the same class of variations so as to have a more like-for-like comparison to the Hjj band.
The Hjj and JetVHeto predictions agree within the Hjj uncertainties, but not quite to within the thickness of the restricted JetVHeto band, in which case the central Hjj prediction is 1-2% below the lower edge of the uncertainty band. Nevertheless, considering the JetVHeto calculation has superior accuracy to both the Hjj and Nnlops predictions, through its high accuracy resummation, the level of agreement we find should be understood as being, again, quite satisfactory: the Hjj prediction is always within 5% of the JetVHeto result, moreover, for the region p T,veto > 25 GeV, it is within 3% of the JetVHeto prediction. We also observed that if we compare to the JetVHeto results with the same uncertainty prescription as ref. [41] (not shown), the central Hjj prediction lies within half the thickness of the more conservative error band that results in that case.
In ref. [29] we presented results showing the Nnlops prediction lying within 1-2% of the JetVHeto prediction, over the full p T,veto range. Some degree of that good agreement stemmed from exploiting freedom in the distribution of the NNLO-to-NLO inclusive K-factor across the leading jet p T spectrum, to 'tune' the Nnlops result. We expect that the slightly less good agreement in the Hjj result here is correlated with the percent level differences seen above in our jet p T spectrum, between Hjj and Nnlops. We remind that these differences are technical in origin, and should be entirely removable with a more refined implementation of our method.
Lastly, we remark that the unimproved Hjj-Minlo results for the jet veto efficiency are, somewhat surprisingly, also quite good. This good agreement of unimproved Hjj-Minlo and JetVHeto is, however, rather fortuitous. The 0-jet cross section in the numerator of the definition of the jet veto efficiency, is equal to the p T integral of the leading jet transverse momentum spectrum from p T = 0 GeV up to p T = p T,veto . One can clearly see from figure 5 that the leading jet transverse momentum spectrum from the unimproved Hjj-Minlo generator is, in general, quite different with respect to the Nnlops and improved Hjj results. For the region p T,veto 30 GeV the Hjj and unimproved Hjj-Minlo jet p T spectra, while clearly different in normalization, are actually not so different in shape. By definition, the jet veto efficiency, ε(p T,veto ), divides out the respective total cross sections, and hence it is therefore reasonable to expect ε(p T,veto ) is not so different in the Hjj and unimproved Hjj-Minlo predictions for the latter p T,veto region. Moreover, since the numerator of ε(p T,veto ) is the cumulant of the leading jet transverse momentum spectrum, which receives, by far, its main contribution from the low p T region, it follows that the behaviour of ε(p T,veto ), for p T,veto 30 GeV, is less sensitive to differences in the latter spectrum in this region, with all predictions converging steadily towards ε(m H ) ≈ 1. We refer back to the discussion of the inclusive 1-jet cross section, surrounding fig. 3, for comments on why one can expect to see small deviations between the Hjj result and the target Nnlops distribution for general inclusive 1-jet quantities, starting at the level ofᾱ 5 S terms. Our initial reaction, to seeing the difference in shape between the y J1 distributions of the Hjj and Nnlops results, was to interpret it as being due to a weakness in our implementation of our method. Re-making this distribution at the level of the Hjj and Nnlops LHE events reveals, however, that the two are actually indistinguishable from one another (the distributions agree at the sub-percent level). Moreover, at the LHE level, the unimproved Hjj-Minlo code is more clearly out of agreement with both of the latter and, in particular, it does no longer agree so well with and the Nnlops in the high rapidity region; the difference being at the level of 5%.
It follows that our implementation of the method actually works perfectly as intended, and that the small features above which were counter to naive expectations, are actually fully attributable to the attachment of the parton shower. The parton shower generates the 3rd hardest radiation and beyond in the Nnlops generator, while it starts by generating the 4th hardest radiation in the case of Hjj-Minlo and Hjj . Naturally then the effect of the parton shower on the y J1 distribution in the Nnlops case is greater, acting to deplete the cross section in the high rapidity side bands relative to the Hjj-Minlo and Hjj results. Given that the difference between the theoretically superior Hjj and Nnlops results in these high rapidity regions has been traced to the effects of the 3rd hardest emitted parton (i.e. anᾱ 5 S effect), we cannot say one result is better than the other. We suffice to say that the difference is in any case small, in a region where theoretical control is not as high as in other places, and it is very much contained within the scale uncertainty bands.  GeV, there is essentially not enough phase space available to generate large L 12 logarithms. By contrast, at the high p J1 T end, one can expect large L 12 logarithms, with a significant contribution from events for which √ y 12 25 GeV and √ y 01 is of the order of p J1 T . So, even though the distribution is defined on 2-jet events, in the high p J1 T limit, the second jet should generally be considered as secondary, soft, radiation emitted from a hard, high-p T , Higgs-plus-jet system. By construction our method will only act to correct the Hjj-Minlo distribution for y 01 y 12 , leaving regions where there is no such strong scale hierarchy untouched. Thus, in fig. 8, at low transverse momentum, we see the Hjj distribution agrees identically with Hjj-Minlo. This is, of course, the desired behaviour, since in this region, for this (2-jet) observable, Hjj-Minlo is nominally NLO accurate, whereas Nnlops is only LO. We remind that, the analogous inclusive leading jet transverse momentum spectrum, fig. 3

, displays significant deviations in shape between
Nnlops/Hjj and Hjj-Minlo in this same p T region, while Nnlops and Hjj are in near perfect agreement.
Turning instead to the high p J1 T region, the three predictions are in good agreement with one another. In the high p J1 T region there is perhaps a faint hint of the Hjj result tending to that of the Nnlops. We assert that the latter tendency would be the correct and desirable result there. Should the transverse momentum of the leading jet enter a high enough p T regime, a 25 GeV jet-defining p T cut for the second jet will correspond to a cut deep in the Sudakov region of the corresponding √ y 12 distribution, in which case, the leading jet p T spectrum in two-jet events increasingly corresponds to the inclusive leading jet p T spectrum.
Finally for fig. 8, we notice that the dark red band, depicting uncertainty due to variations of the technical ρ parameter, has become visible for the first time in this section (in the upper-right ratio plot, at high transverse momentum). This technical systematic is, however, seemingly limited to a ±2% uncertainty, which is dwarfed by the conventional theoretical uncertainty coming from the renormalization and factorization scale variations (the significantly larger light-red band). The last distribution we present showing the behaviour of the leading jet is that of its rapidity in events with at least two-jets, fig. 9. This distribution is rather unremarkable given what we have shown immediately before, for the leading jet transverse momentum spectrum in the same class of events ( fig. 8). Here, as in fig. 8, the distribution shows that the Hjj distribution overlaps the Hjj-Minlo prediction, which is NLO accurate in the descriptions of this observable, while the Nnlops result is only LO. This is the expected and, of course, the desired behaviour of our improved Hjj simulation.

Second jet and third hardest jets
In this subsection we move to present plots of distributions probing directly the behaviour of the second and third hardest jets produced in association with the Higgs boson. As before, jets have been defined according to the anti-k t clustering algorithm, with the jet radius parameter R = 0.4.
Additionally, for the case of jet rapidity distributions, in figures 12 and 13, the jets are required to pass a transverse momentum threshold of 25 GeV.
The transverse momentum spectrum of the second hardest jet is plotted in fig. 10. In all simulations, before (not shown) and after showering, the distribution peaks in the bin at 3 GeV ≤ p J2 T ≤ 6 GeV. Moving upwards from the first bin at p J2 T = 0 GeV the Hjj (red) and Hjj-Minlo (blue) predictions start off with a 20% difference, which smoothly and monotonically diminishes, with the two distributions coalescing at p J2 T ≈ 20 GeV. For higher transverse momenta, the Hjj and Hjj-Minlo histograms become indistinguishable from one another. Meanwhile, in the same region, the Nnlops result starts off with a 15% discrepancy between it and the latter simulations, which rises with the transverse momentum. Nevertheless, the Nnlops prediction is within the margins set by all renormalization and factorization scale uncertainty bands. The behaviour of the Hjj and Hjj-Minlo predictions relative to one another is as intended.
In general, the Hjj-Minlo prediction is NLO accurate in the description of p J2 T , and so it is of course desirable that the Hjj tends to that result in regions where Sudakov logarithms at higher orders are not large, i.e. away from the Sudakov peak. 17 In the vicinity of the peak, large logarithms enter at every order in perturbation theory. In this feasibility study we claim to control these large logarithms nominally at just LL/NLL σ accuracy. The improved Hjj prediction works so as to implement unitarity for the 0-and 1-jet inclusive cross sections by ascribing the mismatch there to missing NNLL σ Sudakov logarithms beyond NLO. The increasing difference of Hjj with respect to Hjj-Minlo in the region p J2 T ≤ 20 GeV, up onto the Sudakov peak, roughly reflects this NNLL σ 'profiling' of the ∼10-12% excess in the Nnlops total inclusive cross section over that of Hjj-Minlo  In figure 11 we plot the transverse momentum of the third jet. In this case there is, coincidentally, good agreement of all predictions in the moderate to high p T domain. This is somewhat fortuitous in the context of the Nnlops simulation, since the third jet in that simulation is generated exclusively in the parton shower approximation, whereas in Hjj and Hjj-Minlo it has a 17 In such regions where it is meaningful to quantify accuracy in the context of just fixed order perturbation theory, we remind that the Nnlops prediction for p J2 T is, by contrast, only LO accurate.
matched matrix element-parton shower description. With a view to validating our ideas, what is more relevant is the observation of the relative behaviour of Hjj and Hjj-Minlo. Here we see, essentially, exactly the same trend as found in the case of p J2 T , specifically, identical agreement for p J2 T 10 − 15 GeV, with a steadily increasing excess of Hjj over Hjj-Minlo as one looks towards zero transverse momentum. These aspects are also fully explained and intended, with the same reasoning as for p J2 T . The only slight difference here is that the third jet being, by definition, softer than the second jet, implies that the excess of Hjj over Hjj-Minlo is confined to a slightly lower region of the p J3 T distribution, than one finds in the p J2 T case.    . 3). As with the p J2 T and p J3 T spectra, for modest values of the transverse momentum, the tendency of Hjj to reproduce Hjj-Minlo here is as intended and desired; the latter being NLO accurate for y J2 and LO accurate for y J3 , in contrast to the LO and parton shower accuracy, respectively, afforded by the Nnlops.

Jet rates
In figures 14 and 15 we present differential jets rates obtained from the exclusive k t -jet clustering algorithm with radius parameter R = 1. Figure 14 shows the log 10 √ y 01 and log 10 √ y 12 jet rate distributions, while figure 15 shows log 10 √ y 23 and log 10 √ y 34 . In the upper plots we display the log 10 √ y01 differential jet rate on the left, while on the right we show the various predictions relative to the central improved Hjj-Minlo (Hjj ), Nnlops and original Hjj-Minlo (Hjj) ones, respectively, in the top, middle and bottom panels. In the lower plots we display the corresponding set of distributions for the log 10 √ y12 differential jet rate. In the making of these plots jets have been clustered according to the kt-jet algorithm, with radius parameter R = 1.
The log 10 √ y 01 distribution in figure 14 is equivalent to a plot of the transverse momentum of the leading jet in the event, defined according to the k t -jet clustering algorithm with R = 1.
It is therefore not surprising to find that the results for this distribution have a markedly similar structure to those for the leading jet transverse momentum spectrum in fig. 5; notwithstanding the fact that in the latter case the jets were defined according to the anti-k t jet algorithm, with radius parameter R = 0.4. We therefore refer the reader back to the discussion surrounding fig. 5, for further explanation regarding the features of the log 10 √ y 01 distribution.
The log 10 √ y 12 in fig. 14 is more interesting, since this distribution is directly affected by our proposed Minlo improvement procedure. One can relatively quickly gain an appreciation for the pattern of the results here by noting that there is some reasonable degree of correspondence to be expected between √ y 12 and p J2 T , based on how √ y 12 is defined; if all the jet clustering algorithm did was initial-state clusterings, they would indeed be exactly the same thing. Despite seeming like an over-simplification, it is nevertheless the case that the relative behaviours of the three predictions here are in very good agreement, quantitatively, with that discussed earlier for the p J2 T spectrum ( fig. 10).

As in p J2
T we see an excess of Nnlops with respect to Hjj and Hjj-Minlo of ∼ 12% in the region 20 GeV √ y 12 100 GeV, with the latter pair of simulations in perfect agreement. For √ y 12 20 GeV, as in the corresponding region of the p J2 T distribution, the Hjj and Hjj-Minlo predictions become increasingly separate, with the former increasing over the latter, manifesting the restorative effect of the correction procedure to recover the inclusive 0-and 1-jet Nnlops cross sections. Even the crossing over of the Nnlops and Hjj distributions appears to occur at exactly the same place in the p J2 T and log 10 √ y 12 distributions ( √ y 12 ∼ 9 GeV). In contrast to the p J2 T distribution, the log 10 √ y 12 plot makes it clearer when, and to what extent, the correction kicksin. One can see that the correction turns on smoothly just before the Sudakov peak, starting at log 10 √ y 12 ≈ 1.25, ( √ y 12 ≈ 18 GeV), leading to a 7% increase in Hjj over Hjj-Minlo on the Sudakov peak, and ranging up to 25% at √ y 12 = 1 GeV.
Lastly, this log 10 √ y 12 distribution shows the first real evidence, so far, of some sensitivity in the Hjj results to the technical ρ parameter. The conservatively estimated systematic uncertainty owing to ρ is depicted by the dark-red band, seen superimposed on the light-red band, in the uppermost ratio plot. This sensitivity to ρ is, however, rather contained at the level of ±10 − 15%, moreover, it is basically negligible above √ y 12 = 3 GeV.
Moving on, in the upper half of fig. 15 we have the log 10 √ y 23 distribution. The correspondence of √ y 12 with p J2 T , which helped to quickly understand the log 10 √ y 12 results above, has an analogon here, namely, that neglecting final-state clusterings by the jet algorithm, √ y 23 becomes equal to p J3 T . This analogy continues to appear to hold remarkably well, for describing the features of log 10 √ y 23 in terms of those found in the p J3 T distribution of fig. 11. The arrangement of the three predictions relative to one another, throughout the log 10 √ y 23 distribution, is very much in direct correspondence with what one can see in the p J3 T distribution. For example, all three predictions even cross at the same point in the log 10 √ y 23 and p J3 T distributions: √ y 23 ≈ 50 GeV in fig. 15 and, correspondingly, p J3 T ≈ 50 GeV in fig. 11. As was noted in comparing the p J2 T and p J3 T distributions beforehand (figs. [10][11], the effect of our corrective procedure in lifting the Hjj distribution above that of its 'parent' Hjj-Minlo simulation, in the region log 10 √ y 12 < 1.25, directly percolates into the same lower reaches of log 10 √ y 23 (and also log 10 √ y 34 ). The extent of this lifting in log 10 √ y 12 and log 10 √ y 23 , is quantitatively compatible with that seen in p J2 T and p J3 T , both in terms of its magnitude and the phase space domain over which it occurs; in particular we note that the separation of the Hjj and Hjj-Minlo distributions starts at a very slightly higher value of y 12 than y 23 , the latter being, by definition, smaller than the former. As with the discussion of the preceding jet rate variables and transverse momentum spectra, the effect of the correction procedure is rather modest and it is limited to a region of phase-space for which all-orders large logarithmic corrections are significant.
In the lower half of figure 15 we show the log 10 √ y 34 distribution. In order to have a non-zero contribution to this observable events must contain at least four partons. So, in the case of Hjj and  Figure 15. In the upper-left plot we show the log 10 √ y23 jet rate, while on the right we show the various predictions again as ratios with respect to one another. In the lower plots we display the corresponding set of distributions for the log 10 √ y34 jet rate. Jets have been constructed using the kt-jet algorithm, with R = 1.
Hjj-Minlo this distribution directly probes, for the first time, radiation which is exclusively due to the parton shower interfacing. The distribution is plainly smooth and exhibits no irregularities that might otherwise signal some problem in that interfacing. The same comments apply here as above, in regards to the lifting of the Hjj distribution with respect to Hjj-Minlo, due to the action of our correction procedure on the y 12 distribution and the associated feed-down from that onto the higher multiplicity differential jet rates.
The penultimate set of differential jet rates we wish to present are given in figure 16. Here we examine the key jet rate of interest to our studies, given its role in the proposed correction procedure, log 10 √ y 12 , but now subject to additional cuts in the √ y 01 jet rate variable. These cuts are intended to bring to the fore events for which there is a hierarchy y 12 y 01 and associated large logarithm L 12 . This aspect is indeed manifested in both log 10 √ y 12 distributions in fig. 16 through the Sudakov peak shifting to higher y 12 values. The Sudakov peak in the inclusive distribution of fig. 14 is centred around log 10 √ y 12 = 1 ( √ y 12 = 10 GeV), moving up to log 10 √ y 12 ≈ 1.5 ( √ y 12 ≈ 30 GeV) on imposing the √ y 01 > 50 GeV cut, as shown in the uppermost plot in fig. 16, and further to log 10 √ y 12 ≈ 1.75 ( √ y 12 ≈ 55 GeV) on imposing the √ y 01 > 200 GeV cut. The shifting of the peak to higher y 12 values is a manifestation of the fact that the cuts imply a proportionate increase in the available phase space for high p T emission of the second pseudoparton.  . 10). In the latter distribution the discrepancy increases with radiation hardness, as it does in fig. 16. Technically, the agreement of Hjj and Hjj-Minlo in this limit is also easy to understand, since in these regions L 12 is not large and the Minlo correction procedure 'switches off', with the Nnlops prediction being categorically inferior to Hjj-Minlo there. Specifically, the h(L 12 ) function (eq. 2.32) tends to zero.
Looking towards the Sudakov peak regions in fig. 16, where the great bulk of the cross section is centred, one expects, by virtue of the fact that our method is to return inclusive 0-and 1-jet Nnlops predictions, the Nnlops and Hjj predictions to agree well there. The results in figure  16 support this simple reasoning quite well.
Turning back to the 1-jet inclusive cross section predictions in fig. 3, with a 50 GeV jet p T threshold, there is a relatively small difference between Hjj-Minlo and Nnlops predictions, this implies that, on average, the δ(Φ BJ ) term in eq. 2.31 is very small in the context of that observable, which is in some reasonable degree of correspondence with the cumulant of the first distribution in for this plot the conclusion is less precise, due to the appearance of large statistical uncertainties.
Taking together the following points, we believe we can now conclude that the uncertainty due to our ρ parameter is in general negligible: i.) we took a rather conservative approach to assessing the uncertainty due to ρ, varying it from 1 to 27, ii.) we constructed observables to isolate and expose potential problems owing to ρ, we found no pathologies, and the latter uncertainty only showed up in the very deep Sudakov region, where theoretical control is very limited.
We conclude the presentation of results on jet rates with the y 12 /y 01 distributions in fig. 17.
The latter quantity is precisely that which our Minlo improvement procedure directly modifies, in order to achieve agreement with the inclusive Φ BJ distribution of the Nnlops (see again sects. 2.6 and 3). The three y 12 /y 01 plots in fig. 17 are also in rough correspondence with those for log 10 √ y 12 in figs. 14 and 16. Indeed, the arrangement of the three predictions relative to one another in fig. 16, for the √ y 01 > 50 GeV and √ y 01 > 200 GeV cuts, is essentially the same as that which one finds for the same √ y 01 cuts applied to the y 12 /y 01 distribution in fig. 17. This correspondence is to expected, if one assumes that the bulk of events making the distributions in both cases (log 10 √ y 12 and y 12 /y 01 ) is dominated by those having √ y 01 close to the cut; this assumption is reasonably fair, given that √ y 01 falls off quickly, like the leading jet p T spectrum. Crudely speaking, this makes the denominator of y 12 /y 01 a constant and, by implication, y 12 /y 01 tends to look the same as scaled a plot of √ y 12 . We also consider that the y 12 /y 01 distribution for √ y 01 > 10 GeV bears a considerable resemblance to that of the inclusive log 10 √ y 12 one in fig. 14. This makes sense on the basis that the √ y 01 cut on the former distribution is loose to the point of being no cut at all. This being the case, we refer back to our discussion on the features of the aforementioned log 10 √ y 12 plots, for explanation of the structures in the y 12 /y 01 ones.

Higgs-jet and dijet correlations
In this subsection we move to check observables more sensitive directional correlations between the Higgs boson and jets in the event. Such variables are routinely encountered in experimental analysis relating to Higgs production via vector boson fusion. In leaving the jet rate variables behind, we return again to define all jets according to the anti-k t clustering algorithm with radius parameter R = 0.4, for all of the remaining numerical results in this work.
We start with the azimuthal separation of the Higgs boson and the leading jet, ∆φ HJ , for events containing at least one jet, in fig. 18. The region ∆φ HJ ≈ π is dominated by configurations consisting of a hard underlying Higgs-plus-one jet configuration, accompanied by additional soft radiations.  In figure 19 we display the azimuthal separation of the two leading jets in the uppermost plot and their rapidity separation in the lower one. In figure 20 we have further plotted the invariant mass of the two leading jets, for events in which they are separated by at least four units of rapidity. All of these distributions demand the presence of at least two jets in the final state. From the analysis of our foregoing results, we understand that for a global jet p T threshold of 25 GeV, we can expect that 2-jet inclusive observables, such as these, are dominated by events with no strong hierarchy of scales y 12 y 01 . Consequently, we expect, and we find, that our corrective procedure has no effect, with the Hjj and Hjj-Minlo results being indistinguishable from one another throughout. This is again our desired behaviour given that the Hjj-Minlo prediction is nominally NLO accurate for these observables, while the Nnlops is similarly just LO. Lastly, we add that the same conclusions hold for the m JJ distribution when the |∆y JJ | > 4 rapidity separation cut is not imposed, in particular, the Hjj and Hjj-Minlo results remain indistinguishable, with the Nnlops continuing to exhibit the same relative discrepancy (albeit within a smaller scale uncertainty band).

Jet binned Higgs boson transverse momentum distribution
The final results we present, in figs. 21-22,    itself with the superior (NNLO) Nnlops result in the low transverse momentum domain. On the other hand, as soon as the Higgs boson has reaches a transverse momentum in excess of that of the jet defining p T threshold, we see that Hjj quickly comes into agreement with Hjj-Minlo. This behaviour is also as intended, since in the latter region, momentum conservation combined with the requirement that there be no resolved jets, dictates that the Higgs boson must be considered as recoiling against multiple hard radiations which are widely separated in angle from one another, all with p T < 30 GeV. The latter class of 'hedgehog' configurations is described more accurately by the higher multiplicity Hjj-Minlo simulation.
Turning to the Higgs transverse momentum in the 1-jet events, we see the results we naively expect in the region p H T > 100 GeV, with Nnlops and Hjj in very good agreement. In the region surrounding the peak of the distribution at p H T ∼ 50 GeV, Hjj continues to agree well with Hjj-Minlo, but not quite as nicely as before. The slight excess of the Hjj prediction over the Nnlops around this peak follows the same explanation as for the similarly sized enhancement of the exclusive 1-jet cross section of the former over the latter, in the discussion surrounding fig. 4.
There we explained that our correction procedure led to an enhanced 1-jet exclusive cross section, by acting to recover the inclusive 1-jet cross section of the Nnlops, while maintaining the 2-jet inclusive cross section of Hjj-Minlo; since the 2-jet inclusive cross section of Hjj-Minlo was low with respect to that of the Nnlops, the Hjj 1-jet exclusive cross section therefore had to be high.
Remarkably, on the other hand, we note that for the lowest bin in the N jets = 1 p H T plot, it is in fact natural and correct that the Hjj distribution is found to be in complete agreement with Hjj-Minlo, for in that region the recoil of the leading jet can no longer be balanced by the Higgs boson, and instead extra radiation must be present to this end.   Lastly, we look to the Higgs boson transverse momentum distributions in the exclusive 2-jet events and inclusive 3-jet events, in the upper and lower plots of fig. 22. For both the exclusive 2-jet and inclusive 3-jet p H T spectra, we see that Hjj agrees perfectly with the Hjj-Minlo generator in the low transverse momentum domain. In the high transverse momentum regions we find that all three predictions agree rather well with one another. In the exclusive 2-jet case at high p H T , we can, however, clearly see that the correction procedure has driven Hjj to reproduce Nnlops rather than Hjj-Minlo. We believe that this too is again the desired result and that the Nnlops prediction is superior to that of Hjj-Minlo in this particular kinematic domain. This assumption is based on the fact that in high p H T N jets = 2 events, the leading jet has a transverse momentum which is bounded from below by approximately half that of the Higgs boson, moreover, it will tend to have a transverse momentum close to that of the Higgs. Thus, the p H T ∼ 200 GeV region of the N jets = 2 events p H T spectrum, will be dominated by events with √ y 01 ∼ 200 GeV. Referring back to the log 10 √ y 12 plot of fig. 16, with the √ y 01 > 200 GeV cut imposed, we can then understand that nearly all such events will come with a second p T > 30 GeV jet 'for free', i.e. the p H T spectrum with N jets = 2, for p H T 200 GeV, becomes essentially the N jets ≥ 1 p H T distribution. Hence, we believe the Nnlops/Hjj prediction to be more accurate than Hjj-Minlo in this case.

Conclusion
In this work we have revisited the Minlo and Minlo frameworks. Our main aim has been to address the issue of how to extend the accuracy of existing Minlo simulations up to that of Minlo . We focused on Minlo simulations of B+2-jet production (Bjj), with B a colourless system, as prototypical 'complex processes', however, our ideas are more widely applicable. For the latter generators, which are NLO accurate in the description of B+2-jet (Bjj) inclusive observables, promotion to Minlo accuracy amounts to the requirement that B+1-jet (Bj) inclusive quantities also be recovered at NLO. We have also considered how to go further in this framework, and obtain (N)NLO accuracy for inclusive B-production observables from Bjj-Minlo.
In existing Minlo simulations the two-fold NLO accuracy is obtained by constructing a Sudakov form factor which returns the relevant inclusive NLO Bmj cross sections, differential in the underlying Born phase space, starting from NLO Bnj cross sections, with m = n − 1. While the form factors are explicitly constructed from high-order resummation ingredients, the accuracy of the resummation in the resulting Minlo simulation is the same as in the initial Minlo one. The net effect of the modifications is to carefully unitarize the inclusive cross section, differential in the Born kinematics. This is very similar to the working of Powheg Sudakov form factors in Nlops matching. 18 The latter contain NLL and even power suppressed terms in the exponent in order to recover NLO accuracy, despite being, in general, just LL accurate.
We started in this work by trying to clarify to what extent the Bjj-Minlo simulations already achieve the aforementioned Minlo accuracy, and to see how to improve them in this direction by better understanding the relevant resummation. We used the Caesar formalism to derive a NNLL σ resummation formula for the 0 → 1-k t -jet rate y 01 and, separately, the 1 → 2 jet rate y 12 ; including leading multiple emission corrections in the exclusive k t -algorithm. The NNLL σ formula reveals existing Bjj-Minlo simulations miss NLL σ terms in their y 12 Sudakov form factor exponents, associated to soft-wide angle gluon emission from the underlying Bj state in the kinematic domain where Caesar is valid, y 01 m 2 H . The Sudakov form factors in existing Bjj-Minlo codes also miss the NNLL σ multiple emission corrections in the resummation formula. With these clarifications one could formally improve Bjj-Minlo codes towards Bjj-Minlo , implementing improved Sudakov form factors to that end.
We derived the fixed order expansion of the NNLL σ Caesar formula and from this we showed how our Minlo procedure applied to the Bj(j) NLO computations returns a matched, resummed, NLO accurate jet resolution spectrum. In doing so we also assumed the presence of unknown N 3 LL σ and NNLL σ terms in our initial fixed order expansion formula; the former owe to the limitations of our initial resummation formula truncated at NLO, while we allowed for the presence of the latter in anticipation of a breakdown of the Caesar framework in considering the region y 01 m 2 H later on.
Upon integration over the Minlo jet resolution spectrum, we determine how the distribution of the Born kinematics differs from that of conventional NLO on account of those unknown terms, which were tracked and contained. We demonstrate how such unwanted terms are removed in the original Minlo approach and, based on that, we introduced an approximately equivalent procedure, which promotes the Minlo simulations to Minlo accuracy without the need of analytic expressions for higher order terms in the Sudakov form factor -terms which are in general unknown. To this end we solve the condition that the missing higher order Sudakov contribution must be such that A more important benefit of the proposed extension is that one can begin to make Minlo simulations, merging two units of multiplicity at NLO, or one NLO and one NNLO, without a merging scale, for complex processes, provided the resummation in the Minlo program entering the correction procedure is NLL σ accurate. The main limitation to be faced in practice will occur for processes in which the dimensionality of the underlying Born phase space is high; in this case the numerical determination of the δ (Φ) term (eq. 2.31) which corrects the Sudakov form factor will become challenging. We should be clear though that, broadly speaking, the main objective in this work has been to demonstrate that one can use precision inclusive cross sections to determine approximate corrections to the Sudakov form factor of Minlo simulations, such that they return the correct (N)NLO Born kinematics; this in turn leads to (N)NLO accuracy for arbitrary infraredsafe observables which nominally receive their leading contributions from parton multiplicities lower than that included in the initial Minlo simulation. The precise way we have done this, discussed in the second part of section 2.6, and our implementation of it, is undoubtedly just one option out of many, and can be simply considered as a practical, working, proof-of-concept at this point.
Even in the worst case, should the dimensionality of the Born kinematics become too much, the method here still has the potential to greatly improve results, in approximating the full δ (Φ) term of eq. 2.31 by a carefully dimensionally reduced version of it.
The loose requirement on the accuracy of the Minlo resummation has an additional useful property: the method can be applied in regions of phase-space where the underlying Born itself has disparate kinematic scales associated with it; in such regions achieving high accuracy resummation is currently a formidable challenge. In particular, our improvement procedure remains valid for Bjj-Minlo in regions of phase-space where y 01 is much smaller than m 2 B , where both large logarithms of m 2 B /y 01 and y 01 /y 12 require resummation. To this end we first argued that the Caesar y 01 /y 12 resummation remains NLL σ accurate in the region y 01 m 2 B . Our argument is largely based on the finding that the Caesar Sudakov form factor for this variable is the same as that prescribed by the coherent branching formalism at NLL σ (except for a subset of the soft-wide-angle radiation, which is beyond the accuracy of the latter). We also note that the leading NLL σ terms in the fixed order expansion of the Caesar resummation formula are obtained by integrating a single soft/collinear emission over the region y 12 < y 01 , i.e. they are the same whether y 01 m 2 B , or not. Taking the latter two points together, it follows that the Caesar y 01 /y 12 resummation must hold at the NLL σ level, even in the difficult regions.
The fact that y 12 resummation works in all regions of phase-space implies one can do a 'nested' Minlo simulation with the help of our proposed extension: instead of using unitarity such that partially integrated Bjj distributions become equal to the NLO Bj distributions, we train them on Minlo Bj distributions. This makes them also NLO correct in inclusive B distributions, since the latter are obtained on suitably integrating over radiation in Bj-Minlo . In this way the extension of the Minlo method to the merging of more than two multiplicities is realised.
As a feasibility study for the latter, we have applied our correction procedure to Hjj-Minlo.
We start from a LHE event file for the Hjj-Minlo simulation and reweight the events such that distributions differential in Φ HJ become equal to the existing NNLO-improved Hj-Minlo calculation, without hampering the formal accuracy of the Hjj-Minlo simulation. We therefore made predictions that are NNLO accurate for inclusive Higgs boson production, and NLO accurate for Hj and Hjj observables. Since the LHE events are ultimately generated according to the Powheg Nlops matching procedure they may, of course, be showered in the usual way. Our numerical results are very encouraging: for inclusive observables in H and Hj production, we recover the results of the Hnnlops simulation, while for observables in which y 12 ∼ y 01 we recover the Hjj-Minlo predictions, with smooth interpolation between them.
There is ample freedom in the functional form of the reweighting factor which is formally beyond the accuracy of the method. We have explored (some of) its dependence and seen essentially no visible effects of it in the many distributions we have examined. The distribution which displayed most sensitivity to this ambiguity was, unsurprisingly, log 10 √ y 12 . Even for this variable, the sensitivity is located in the deep Sudakov region, mostly well-below the Sudakov peak, a region which is anyway very sensitive to higher order resummation and non-perturbative effects.
There are a number of aspects of this work which can be explored further and refined. It is clear, for example, that it is interesting to consider our approach in application to other processes.
We have shown the method can work well for a process with 3 final-state particles (Hjj), thus it seems reasonable to expect similar quality results in application to processes with equal multiplicity, e.g. trijet, and jet-associated single-top/top-pair production (with some approximation in the handling of top decays). In fact these processes are in one sense less demanding than that which we demonstrated, in so far as we dealt with a process for which two jets could become unresolved, moreover, this was handled while mapping onto an NNLO calculation of Higgs production. On the other hand, for high multiplicity processes like VBF Higgs-plus-3-jet production, the dimensionality of the phase space combined with the problems to be anticipated in obtaining high statistics for determining δ (Φ), would likely prove too cumbersome in practice, at least for our proof-of-concept implementation. Nevertheless, in the absence of a better alternative, we would still advocate trying the latter method in some approximate form, e.g. applying it on only a carefully selected subset of the variables which parametrize the underlying Born kinematics. Depending on the initial circumstances this may lead to very desirable improvements.
Relatedly, on a technical level it is worth considering a more sophisticated approach to our im-plementation, e.g. using an adaptive, optimised, grid parametrization procedure for the underlying Born phase space, together with advanced interpolation procedures for computing δ (Φ). Having said that, it is perhaps a good indicator of the potential of this approach, that it appears to have worked remarkably well even with just a basic implementation using hand-made rigid grids. From a theoretical perspective, one may wish to consider improvements to the Bjj-Minlo codes based on our comparison of their inherent resummation and that of the Caesar framework, such as inclusion of soft-wide-angle (∆S 1 ) and multiple emission (F 2 ) terms in the Sudakov form factors.
Our numerical studies in this paper suggest that these inclusions would be of really quite limited interest though.
A further investigation would be to consider the effect of breaking the reweighting procedure for Bjj-Minlo into two phases: in the first stage just the inclusive Φ B distribution of Bjj-Minlo is corrected to that of Bj-Minlo , by adjusting the y 01 distribution; in the second stage the procedure is applied to the Bjj-Minlo output from the first stage, in exactly the same way as set out in sects. 3-4. At the NLL σ level, there is no distinguishing between the latter approach and that which we carried out. On the other hand, it is clear that, at some level, the effective Sudakov form factor correction that we derive for the y 01 /y 12 resummation will make up for what might better be considered as deficiencies in the y 01 Sudakov form factor. Nevertheless, from our numerical studies here, we expect that this change would only register much like the ρ-parameter variations that we assessed, i.e. we believe it will only become visible in y 12 regions which are under poor theoretical control (the deep Sudakov region). It is also important not to over emphasise this point in view of the fact that the correction procedure obtains the Φ B , y 01 and y J distributions of Bj-Minlo , by construction, in any case. Nevertheless, this alternative may prove advantageous in other applications of the method.

A Appendix
A.1 Derivation of LL/NNLL σ jet resolution spectra In this appendix we give details on how our general NNLL σ jet resolution spectrum is arrived at from the results of refs. [55] and [59].
In eqs. 3.14-3.17 of ref. [59] the NLL resummed cumulant cross section, matched to NLO, to yield also NNLL σ accuracy is given as with dσ 0 /dΦ the leading order cross section fully differential in the Born kinematics (denoted B in [55,59]), and with f (v) encoding the resummation. The term α S C 1 is the matching coefficient defined by (eq. 3.16 of ref. [59]) where the dΣ NLO (L) in our notation corresponds to dΣ 1 (v) of [59], and with our dΣ R,1 (L) corresponding to dΣ r,1 (v) of [59]. Thus the matching coefficient of [59] is in our notation α S C 1 = α Sχ1 (Φ). The function f (v) is the main result of ref. [55] (see eq. 3.6 therein) and is comprised as follows (taking b = 0, a = 2, d = g = 1, and henced = 1, as appropriate for the k t -jet resolution variables considered in our work, namely, V ({p} , k) = (k into that for f (v) gives, with no approximations, with −R (v) here as given in our eq. 2.4. Now we start to make approximations and deviate from ref. [55], breaking the NLL resummation by N 3 LL σ terms. Consider that F (R ) resums single log terms as F (R ) = 1 + F 2 R 2 + ... + O (F n R n ) + ... , (A.7) R = ∂ L R (v) and so R n is O (ᾱ n S L n ). Neglecting terms of N 3 LL σ accuracy we can simply replace R =ᾱ S (y) 2G 12 L F (R ) = exp and hence We then have the LL/NNLL σ resummed cumulant expression . (A.10) Recall thatχ 1 encodes hard-virtual and hard-collinear splitting corrections, and that these contributions contain terms which cancel the µ R and µ F dependence of dσ 0 /dΦ B . We may separate these parts as follows:χ where H 1 Φ, µ 2 R , Q 2 = H 1 (Φ) + 2qβ 0 ln µ R Q + G 11 + 2S 1 − 2G 12 ln Q qq Q − qβ 0 2 ln Q qq Q ,(A.12) and C 1,ij z, µ 2 F , Q 2 = C 1,ij (z) − 2 ln µ F Q P ij (z) , C 1,ij (z) = −P ij (z) − δ ij δ (1 − z) A ij π 2 12 . (A. 13) We underline that in the relation betweenχ, H and C 1 , eq. A.11, the renormalization scale in H 1 is set to µ R and in C 1 , which is convoluted with a PDF, the explicit factorization scale therein is µ F , i.e. C 1 in eq. A.11 is precisely as it is written in eq. A.13. Equation A.12 basically defines H 1 as what remains of H 1 when its Q and µ R dependence is removed, q being the number of powers ofᾱ S associated to the Born process. The P ij (z) terms are the O ( ) parts of the LO splitting function P ij (z) : P qq (z) = −C F (1 − z) , A qq = C F , P gq (z) = −C F z , A gq = 0 , P qg (z) = −z (1 − z) , A qg = 0 , P gg (z) = 0 , A gg = C A . (A.14) Equations A.11 and A.12 serve to define H 1 . Whereas H 1 is dependent on the virtual corrections to the underlying hard scattering process, the C 1 pieces are due to collinear splitting and only depend on the flavour of the (incoming) legs in the Born configuration. While we strictly only aim for NNLL σ accuracy in our initial resummation formula, to better enable comparison/extension with existing NNLL work, without affecting any NNLL σ terms, we opted to replace in our resummation formula (eq. A.10) .

(A.15)
Note in particular the introduction of the v dependence in the renormalisation and factorisation scales in the final term. From here eq. 2.1 follows immediately on differentiation with respect to L.

A.2 Fixed order expansion of resummation and Minlo formulae
Here in eqs. A. 16 The factorization scale in C 1 in eqs. A.18-A.19 (including the C 1 implicit inχ) is set to µ F ; it is exactly as written in eq. A.13. We point out that for the regime y 01 m 2 B , in Bjj-Minlo, the virtual corrections, H 1 , will contain large logarithms of ratios of scales deriving from the related underlying (Bj) Born kinematics approaching a singular region.
The dσ S expansion of eq. 2.20, with H nm as given in eqs. A. 16-A.19, is invariant under µ R and µ F variations up to higher order terms (∝ dσ 0ᾱ 3 S ) beyond NLO accuracy. Also dσ S with these H nm is invariant under variations of the resummation scale, Q, up to and including NNLL σ terms. To make dσ S invariant under resummation scale variations also at the N 3 LL σ level, requires modifying H 20 → H 20 + [K + 4F 2 G 12 ] H 11 ln(Q 2 qq /Q 2 ) only. Such a modification could be easily generated by simple adjustment of our initial resummation formula, however, since the latter is only guaranteed to reproduce terms to NNLL σ accuracy anyway, we do not consider this.
We have compared our expansion formula, eq. 2.20, to known results for the W/Z and Higgs boson transverse momentum spectra [35,37], as well as to those of the associated leading jet (derived by expanding the NNLL resummation of Banfi et al [41]). To ease comparisons, we note the following relations between our notation and refs. [35,37,41]: where A (1) , B (1) and A (2) are used in the latter articles. We also point out that for B production F 2 = 0 in both the B transverse momentum spectrum and that of the leading jet (eq. 2.10). Lastly, in the results of refs. [35,37] the resummation scale is set to the invariant mass of B, i.e. Q = Q qq in our notation, leading to S 1 → 0 here, as well as simplifications in theχ and H 1 functions.
With the correspondence in notation understood, we find complete agreement between our singular NLO expansion formula, eqs. 2.19-2.20, and those of refs. [35,37,41], up to and including NNLL σ terms. To also have agreement with refs. [35,37] for the N 3 LL σ terms in the Higgs/vector boson transverse momentum spectrum, we need only add to H 20 , in eq. A.19, with B (2) as given in ref. [41]. For full agreement with ref. [41], including N 3 LL σ terms, we only have to add to H 20 , in eq. A.19, where f clust and f correl are corrections due to jet radius dependent clustering/correlated emission effects, as given in ref. [41]. We point out that the needed/missing R 20 term here (eqs. A.21-A.22), for B-production, is just a number with no dependence on kinematics.
(A. 25) 4. Multiply by the Sudakov form factor timesᾱ s K 2 R y /ᾱ s µ 2 R : Precisely, the steps outlined above are those used in the construction of the Bj-Minlo simulation of ref. [28], adopting the general notation of sect. 2, so that they apply equally to Bjj-Minlo We are interested in the expansion of the latter up to and including NLO terms. Noting that dσ R is a total derivative with respect to L, and using the definitions ofχ in terms of H 1 and C 1 , we have dΣ M,1 (L) dΦ = dΣ R,1 (L) dΦ In determining the equality in eq. A.28 we have made use of the relation in eq. 2.14. In going from eq. A.28 to eq. A.29 we have made use of theχ 1 definition in eq. 2.15.

A.5 δ (Φ) denominator
Neglecting N 3 LL σ terms (as in sect. 2.6) we obtain for the denominator of δ (Φ) in eq. 2.31 where the O ( √ᾱ S ) comes from N 3 LL σ terms in the integrand.