Resummation ambiguities in the Higgs transverse-momentum spectrum in the Standard Model and beyond

We study the prediction for the Higgs transverse momentum distribution in gluon fusion and focus on the problem of matching fixed- and all-order perturbative results. The main sources of matching ambiguities on this distribution are investigated by means of a twofold comparison. On the one hand, we present a detailed qualitative and quantitative comparison of two recently introduced algorithms for determining the matching scale. On the other hand, we apply the results of both methods to three widely used approaches for the resummation of logarithmically enhanced contributions at small transverse momenta: the MC@NLO and POWHEG Monte Carlo approaches, and analytic resummation. While the three sets of results are largely compatible in the low-pT region, they exhibit sizable differences at large pT. We show that these differences can be significantly reduced by suitable modifications of formally subleading terms in the Monte Carlo implementations. We apply our study to the Standard Model Higgs boson and to the neutral Higgs bosons of the Two-Higgs-Doublet Model for representative scenarios of the parameter space, where the top- and bottom-quark diagrams enter the cross section at different strength.


Introduction
After the discovery of a Higgs boson at the LHC in 2012 [3,4], many detailed studies have analyzed its properties in order to assess its compatibility with the Standard Model (SM) (see, e.g., refs. [5,6]). These analyses rely on accurate theoretical predictions of Higgs cross sections and decay rates. Recently, the first experimental results for the Higgs transversemomentum (p ⊥ ) distribution have been published in ref. [7,8]. Such measurements, in particular with the increased statistics expected from LHC Run II, open up many new interesting possibilities to test the nature of the Higgs couplings to the SM fields, and thus to probe for signs of physics beyond the SM (BSM). Indeed, the gluon-gluon-Higgs vertex could be mediated by loops of non-SM particles which can affect the shape of the -1 -

JHEP01(2016)090
denote the Higgs and the top-quark mass, respectively), so it seems obvious that for the dominant top-loop induced contribution the matching scale should be chosen to be of the order of m h . However, despite its rather small contribution to the total cross section, the bottom-loop plays a non-negligible role at small values of the Higgs p ⊥ . Since it involves two very different scales (m h and the bottom-quark mass m b ), the choice of the matching scale is not at all obvious in this case. In fact, it was observed that the theoretical prediction depends quite sensitively on the resummation method as well as on the nature of the respective matching scale [11,48,49].
Recently, two algorithmic strategies for determining an adequate central value for the matching scale have been proposed [1,2]. It is important to note that, even though the matching scales of the various approaches have different origins and meanings, i.e., they deal with the unknown higher-order terms in different ways, the proposed strategies may be applied to all of the resummation and matching approaches mentioned above, as will be discussed in more detail below.
From this discussion, it is clear that the resummed transverse-momentum distribution suffers from a number of ambiguities which are formally of higher logarithmic order. It is the goal of this paper to study their numerical impact on the Higgs p ⊥ distribution in and beyond the SM. Since the explicit analytic form of these ambiguities is not fully accessible, our analysis will rely on the numerical comparison of the results obtained in the various approaches, assuming specific phenomenological parameters.
The implementations of these approaches on which we base our study are all publicly available: • for the AR approach, we use MoRe-SusHi, which includes the description of the resummed p ⊥ distribution at next-to-leading logarithmic (NLL) accuracy consistently matched to the fixed-order cross section at NLO QCD [1,48,51]; • the NLO+PS accurate POWHEG implementations of the gluon-fusion process are contained in the directories gg H quark-mass-effects and gg H 2HDM [11] of the POWHEG-BOX [52,53]; • the corresponding MC@NLO implementation at NLO+PS is available in aMCSusHi [54,55] which combines the MadGraph5 aMC@NLO package [56] with the SusHi amplitudes [57].

JHEP01(2016)090
Leaving aside the specific values of the associated matching scales, we will refer to these three approaches and their implementations as "resummation codes" or simply "codes". All codes work at NLO QCD accuracy in the prediction of the Higgs production total cross section, i.e., O(α 3 s ). The differences in the p ⊥ distribution are formally subleading, 5 but can be numerically sizable, as we will see later on. In order to assess the impact of these differences, we compare their numerical results using the same values of the matching scales for all of them. On the other hand, we compare the results of a single code for the two different strategies of setting the matching scale proposed in refs. [1,2]. As we will see, both the intrinsic difference in the formulation of the codes as well as the dependence on their matching scales are a source of sizable ambiguities in the theoretical prediction of the Higgs p ⊥ distribution, in particular at intermediate and large p ⊥ . The sources of these differences will be investigated in detail in the course of this paper.
In the SM, the matching of fixed-order and resummed results has been achieved with NNLO QCD accuracy on the top-quark induced component of the total cross section, i.e., NNLO+next-to-NLL (NNLL) in AR [40,49,58,59] and NNLO+PS accurate Monte Carlo [60][61][62][63][64]. Since our main focus is on the p ⊥ distribution in BSM scenarios with large bottomquark effects and the possibility of new additional heavy states, we refrain from including such effects in our discussion.
This paper is organized as follows: in section 2 we introduce the gluon-fusion process at NLO QCD and summarize the underlying resummation procedures of the three codes under consideration. Our focus is on the respective matching prescriptions and corresponding matching scales. In section 3 we recall the two strategies of refs. [1,2] to determine proper values for the matching scales. The two approaches are then subject to a qualitative and quantitative comparison, where we quote values for the matching scales in a large range of Higgs masses. Our main study of the matching ambiguities is presented in section 4 for the SM and for various scenarios in a generic type-II Two-Higgs-Doublet Model (2HDM), designed in order to emphasize specific contributions to the gluon-fusion Higgs cross section. Section 5 contains our conclusions.

Resummation procedures
In this section, we briefly review the ingredients that are required for the theoretical prediction of the Higgs production cross section via gluon fusion with NLO accuracy. We then summarize the main features of three procedures that allow to resum logarithmic terms in the Higgs transverse-momentum distribution at small p ⊥ , namely AR, MC@NLO, POWHEG. The main focus will be put on their matching prescription and the corresponding matching scales.
Concerning the total Higgs production cross section, the three codes considered in this paper work at the same perturbative accuracy, NLO QCD, using the same matrix elements for the description of the virtual corrections and the real-emission effects. We introduce a 5 Note that the meaning of "subleading terms" is somewhat different for AR and the MC generators. AR consistently resums NLL terms to all orders, while the PS in the Monte Carlo approaches strictly includes only the leading logarithms, but resums also some logarithms beyond the leading ones. common notation to identify these contributions in the three different approaches for the matching of resummed and fixed-order results. The Born-level squared matrix elements B determine the LO cross section and multiply all the soft and collinear counter terms. The corresponding loop-induced Feynman diagrams are shown in figure 1 (a), with both the top and the bottom quark running in the fermionic loop. 6 The interference of the UV-and IR-regularized virtual corrections with the Born amplitude will be denoted byV fin in the following. 7 Its evaluation requires the computation of the two-loop virtual diagrams [65][66][67][68], e.g, in figure 1 (b). Note that there are some differences in the definition ofV fin , due to the IR-regularization adopted in each matching procedure. We simply assumeV fin to be properly subtracted in the respective approach under consideration. At the same perturbative order, real-emission subprocesses need to be taken into account [65,69]. They involve an additional final-state parton with respect to the Born-level process; their squared matrix elements are collectively called R. Some examples of the corresponding diagrams are shown in figure 1 (c-d).
For convenience, we use the following symbolic notation for the convolution over the parton-density functions (PDFs): where M ij (q 2 , µ 2 F , z 1 , z 2 ) is the squared matrix element for the scattering of partons i and j carrying proton momentum fractions z 1 and z 2 , respectively, q 2 is the momentum transfer of the scattering process, is the product of two PDFs, and µ F is the factorization scale.

Analytic p ⊥ resummation (AR)
The first resummation procedure that we consider is the analytic resummation of soft and collinear logarithms in the inclusive transverse-momentum spectrum, see refs. [38,40].
For the matching of the low-and high-p ⊥ Higgs cross section at NLO+NLL, 8 the additive procedure of ref. [40] is adopted, and the hadronic differential cross section is 6 Diagrams courtesy of S. Liebler. 7 We use a hat to indicate IR-regularized quantities in this paper. 8 Note that the accuracy (including terms of order α 3 s ) at large p ⊥ is formally only LO.

JHEP01(2016)090
obtained as dσ dp 2 where dΦ B and dΦ represent the phase space of the Born process and of the process that includes one additional real parton, respectively. By dX/dp 2 ⊥ we denote integration over all variables of X except p 2 ⊥ . The function F is defined as follows: is a numerical constant, 9 J 0 (x) is the Bessel function of the first kind with J 0 (0) = 1, S is the hadronic center of mass energy, and the sum over i, j runs over all kinds of partons.
In the notation of the present paper, we apply the "hard scheme" as defined in ref. [70] for the collinear coefficient functions C ab (z); their first order expressions for gluon-induced processes can be also found in that reference. The Sudakov form factor S(α s ,L) accounts for the resummation of logarithms of the formL = ln(b 2 Q 2 res /b 2 0 +1) at NLL accuracy, with α sL being considered of order unity. The functions g (1) and g (2) , relevant for leading logarithmic (LL) and NLL accuracy, respectively, are given in ref. [40]. The scale Q res is conventionally called resummation scale and controls up to which values of p ⊥ the resummation is effective. Since it thus parameterizes the arbitrariness in the separation of the "soft" from the "hard" region, it plays the role of the matching scale in the AR framework and is usually chosen at the order of characteristic scale of the hard scattering process. Note that the first and the third term on the r.h.s. of eq. (2.2) are explicitly Q res dependent. However, this dependence is canceled through the matching formula order by order in the logarithmic expansion, which renders the p ⊥ distribution independent of Q res when computed to infinite logarithmic accuracy. In eq. (2.2), F NLO is the NLO truncation of F NLL ; i.e., for gluon fusion, it includes terms up to α 3 s . The last term in eq. (2.2) subtracts the singular behavior of the real corrections from the fixed-order expression, given by the second term on the r.h.s.; it therefore avoids double counting of logarithmic terms which are already contained in the first term.
Furthermore, the matching procedure of eq. (2.2) induces unitarity on the matched cross section, meaning that integration of eq. (2.2) over p 2 ⊥ yields the total NLO cross section: dp 2 ⊥ dσ dp 2 The AR approach has been implemented, using parts of HqT [40,58,59], in the code MoRe-SusHi [1,48,51] that will be used in all the numerical simulations based on AR in this paper.

NLO+PS Monte Carlo
An alternative to the analytic resummation is offered by PS Monte Carlo event generators, where PS algorithms allow the numerical simulation of multiple parton emissions. A consistent matching of the fixed-order NLO QCD predictions 10 with the PS has been discussed in ref. [46] and in refs. [47,50,52], and implemented in the MC@NLO Monte Carlo event generator [71] and in the POWHEG-BOX [72], respectively. We refer the reader to the above publications for a detailed discussion of the two implementations and summarize here the main differences of these two approaches, from the point of view of the matching of resummed and fixed-order results.
We introduce the "shower operator" I n (t 1 ) to represent in a compact form the action of the PS to describe n parton emissions; the latter are ordered with respect to a parameter t, which for simplicity we assume to be the transverse momentum of the emission in this paper (other options are the virtuality, or the angle, for example), starting from a maximum value t = t 1 . For the i th emission of a parton, the PS associates a Sudakov form factor times an approximate emission probability, both evaluated at the same value t i of the specific ordering parameter. The ordering of the emissions is a requirement for the PS algorithm to reach the LL accuracy.
In both the MC@NLO and POWHEG approaches, the hardest parton emission is treated retaining the accuracy of the exact matrix elements, whereas the others are generated according to the PS algorithm.
Given a hard scattering process, we describe the evaluation of the cross section dσ n for the radiation of n additional partons in terms of two steps. The first step results in the cross section dσ 1 , which includes only the hardest emission; the remaining n − 1 emissions are taken into account in the second step. This splitting is represented symbolically as applying I n−1 (t 1 ) to dσ 1 , and multiplying by the relevant phase space of the additional n−1 particles: dσ n = I n−1 (t 1 ) dσ 1 dΦ n−1 . (2.5) The formula that describes the hard scattering with the emission of 0 or 1 additional parton can be written, in a sufficiently general way, as 11 where t min is the value of the ordering variable t below which no emissions are allowed, 12 and where we used the factorization of the real phase space into the Born one times the one of a single additional real parton, dΦ = dΦ B dΦ r . 10 Note also in this case that the description is formally only LO accurate as far as large p ⊥ are concerned. 11 While in the POWHEG approach the first emission is explicitly implemented as given in eq. (2.6), it is vital for the subsequent discussion to stress that, in the case of MC@NLO, the terms in the curly brackets in the first term of eq. (2.6) are implicitly generated when attaching the shower to the Born-level configurations. 12 The physics below the scale tmin is not properly accounted for by a perturbative description. Its numerical value is typically set to the hadronization scale of the PS, coherent with the fact that the first term in the curly bracket describes the no-emission probability.

JHEP01(2016)090
As defined before, R contains collectively all the squared matrix elements of the realemission subprocesses, and we further split them into two groups: R div is the sum of squared matrix elements which are divergent in the limit of soft or collinear emissions (in our case gg → gφ and qg → qφ); the regular ones are denoted by R reg instead (in our case qq → gφ). The squared matrix elements of the divergent subprocesses can be further split in two parts: The term R s contains the soft-and collinear-singular part of R div , while R f is the finite remainder. 13 Obviously, this splitting is not unique, since finite parts can be shifted between R s and R f . The generalized Sudakov form factor [73] is denoted by the symbol ∆ s t in eq. (2.6), with t the shower ordering variable; it depends on R s and expresses the probability of not emitting any parton with a value for the ordering variable larger than its own argument t: The factorB in eq. (2.6) is defined bȳ and includes the contributions of the Born squared matrix elements, the corresponding virtual corrections, and the integral over the radiation phase space of R s . The finiteness of theB factor is guaranteed by the fact that all the divergent terms are properly subtracted; this is possible thanks to the renormalization of the UV divergences, to the cancellation of the IR soft singularities between virtual and soft real contributions, and to the cancellation of the collinear singularities, reabsorbed in the definition of the physical proton PDFs. The curly bracket in eq. (2.6) describes the probability of zero or one parton emission in those subprocesses that are divergent in the soft/collinear limit, where R s is the singular part of the squared matrix elements. The precise definition of R s (or, equivalently, R f ) therefore directly affects the expression of the Sudakov form factor. In the following we comment on the two choices adopted in the POWHEG and in the MC@NLO implementations. We stress that these two alternatives differ by terms that are formally of higher order in the perturbative expansion, but that can nevertheless be numerically sizable. The arbitrariness in the definition of R s can be exploited to parameterize the matching procedure uncertainties.
The last two terms in eq. (2.6) depend on the process (R reg ) and on the definition adopted to split R s,f ; both yield a regular contribution in the soft and the collinear limit.
The evaluation of the exact real and virtual matrix elements guarantees, for any observable inclusive over radiation, the NLO QCD accuracy. The latter is preserved by the unitarity of the PS algorithm in the generation of each additional real parton; this feature holds also for the first emission, described by the curly bracket in eq. (2.6). 13 The splitting into R div and Rreg seems redundant at this point, but will play a role in the POWHEG approach, see section 2.2.2.

MC@NLO
In the MC@NLO formulation, the PS algorithm is used to generate all the additional parton emissions starting from the Born-level and real-emission configurations. The exact O(α s ) corrections are applied in order to recover the exact matrix element description of the first hard emission and the correct normalization, including the effect of the virtual corrections to the underlying Born.
The Sudakov form factor implemented in the PS generators uses a universal, processindependent expression to describe parton radiation in the soft and collinear limit, based on the Altarelli-Parisi splitting functions P . Using the notation of eq. (2.7), in MC@NLO the singular part of the squared real-emission matrix element is R s MC@NLO ∝ α s P B, and the generalized Sudakov form factor that appears in eq. (2.6) is actually not explicitly implemented, but generated by the first PS emission on top of the Born-level configuration.
Given the assignment for R s MC@NLO , the definition of R f MC@NLO = R div − R s MC@NLO follows by construction, as the difference between the exact real matrix element correction and its PS approximation, sometimes called Monte Carlo subtraction term; the latter is needed to avoid a double counting with the emission described by the second term in the curly brackets of eq. (2.6) that is generated by the first emission of the shower.
In the MC@NLO approach, the differential NLO+PS cross section with respect to a variable O is: where the sum runs over all possible n-parton configuration after the shower in the final state. We shall stress at this point that the shower spectrum I n in the first and I n−1 in the second line start from different configurations, i.e., from Born-level and real configurations, respectively, and that the observable O is defined on a different phase space in the two lines. The superscript MC is attached when the phase-space is computed in the PS approximation. The real Monte Carlo phase space dΦ MC tends to dΦ in the IR limits, and dΦ MC r = dΦ MC /dΦ B is therefore the PS approximation of the one-particle phase space. By dX dO we denote integration over all variables of X except O.
In the first line of eq. (2.10) all emissions are described via the PS, denoted as "soft" events, while in the second line we have the so called "hard" events: the exact matrix elements with one additional real emission with respect to the Born are sampled over the real phase space; to avoid a double counting with the first line, there is a MC subtraction term, which is evaluated over the approximated one-particle phase space and renders the expression in the squared brackets finite. Expanding eq. (2.10) up to the first emission (order α s with respect to the Born), we indeed recover the structure of dσ 1 in eq. (2.6).
In the MC@NLO approach, the shower emissions are bounded from above by identifying t 1 ≡ Q sh . The so-called "shower scale" Q sh therefore determines the hardest scale accessible to the shower emissions and is usually set to the characteristic scale of the hard scattering -9 -JHEP01(2016)090 process. Q sh plays the role of the matching scale, since it separates the soft/collinear from the hard region in the additive MC@NLO matching approach -indeed very similar to the role of the scale Q res in AR. More precisely, one may choose different shower scales for the soft (Q s sh ) and the hard (Q h sh ) events, as indicated in eq. (2.10). A general feature of the default MC@NLO formulation is that the strict bound enforced by Q sh on the shower emissions suppresses the PS contribution at values of p ⊥ Q sh , resulting in the recovery of the fixed-order distribution at sufficiently large p ⊥ , which is regarded to provide a proper prediction in that region.
The original choice for Q sh in the MC@NLO code [46,71] was to set its value to the shower scale of the PS Monte Carlo, which corresponds to a narrow distribution around the Higgs boson mass. In MadGraph5 aMC@NLO, on the other hand, the shower scale for the soft events, Q s sh , is statistically extracted from a probability distribution with a range that can be defined by the user. 14 The shower scale for the hard events, Q h sh , is fixed to a specific value (i.e., using a δ-distribution) that by default is chosen to be the upper bound of the probability distribution applied to the soft events. For reference, figure 2 shows an example of the shower scale distribution of all events for a SM Higgs boson (m h = 125 GeV) with the restriction Q sh /GeV ∈ [11.425, 114.25]. The spike at the upper end of the distribution is due to the hard events. The peak at the center of the distribution can be considered as "effective" shower scale; we will often refer to it simply as "shower scale" in what follows. E.g., figure 2 has an (effective) shower scale of Q sh = m h /2 = 62.5 GeV. In this paper, we always use a distribution with a shape as shown in figure 2, centered around the respective matching scale, and with ratio of the endpoints equal to ten [54]. The MC@NLO results in this paper have been obtained with the code aMCSusHi [54,55] which combines the MadGraph5 aMC@NLO framework with the matrix elements of the code SusHi.
14 Further details can be found in ref. [56]; this technique has been used sparingly also in MC@NLO v3.3 [74] and higher.

POWHEG
The fully differential cross section of an event with the emission of additional partons is obtained by inserting eq. (2.6) in eq. (2.5). For a generic observable O we have where the shower scale t 1 is set equal to the transverse momentum p rad ⊥ of the parton radiated using the POWHEG approach. 15 In the POWHEG formulation, the splitting of eq. (2.7) is obtained through a dynamical, i.e. p ⊥ -dependent, damping factor D h (2.12) In this case, the role of the matching scale is assumed by h. For Higgs transverse momenta larger than the scale h, the damping factor suppresses R s , while the Sudakov form factor quickly approaches 1 and the spectrum is described by the finite remnant R f POWHEG . Instead, when p ⊥ → 0, R s POWHEG tends to R div . In this limit, R div factorizes into the product of the underlying Born multiplied by the universal Altarelli-Parisi splitting functions, and the POWHEG Sudakov form factor yields a suppression of the transverse-momentum distribution. The role of the scale h can be understood from two points of view: i) it is the maximum value of Higgs transverse momenta for which the curly bracket in eq. (2.6) is appreciably different from zero; the normalization factorB multiplies the curly brackets and rescales them in this p ⊥ interval; ii) considering that the (POWHEG) Sudakov form factor is a function that varies between zero and one, the scale h controls the region of p ⊥ where the suppression is active.
A general feature of the POWHEG approach is that it generates a tail of the p ⊥ distribution that is higher than the fixed-order prediction in this region. The description of the enhanced p ⊥ tail and, in particular, of the weight assigned to each high-p ⊥ event, deserves some discussion. The emission probability given by the squared real-emission matrix elements R in eq. (2.6) is proportional to α s and is multiplied by the overall factorB, which starts at LO but includes also O(α s ) corrections. The latter are related to the total K-factor and enhance the first emission weight. 16 The damping factor D h allows to reduce the p ⊥ interval where this reweighting is active. A second relevant element to understand the enhancement of the large-p ⊥ tail is given by the impact of multiple emissions beyond the first one. In POWHEG the phase space available for the second emission is limited only by the p ⊥ value of the first one and changes on an event-by-event basis. In some cases, the second and following emissions can still be very hard. This explains why the POWHEG high-p ⊥ tail tends to be JHEP01(2016)090 larger than the one obtained in the other two approaches and at fixed order. This formulation differs in particular from the one adopted in MC@NLO, where instead all the PS emissions are limited by the shower scale, so that the p ⊥ distribution merges into the fixed-order one at sufficiently large transverse momenta. We will return to this discussion in section 4.2.3.
The POWHEG results in this paper have been obtained with the gg H quark-mass-effects and gg H 2HDM generators [11]. They are both part of the POWHEG-BOX framework [53].

Determination of the matching scale
In this section we describe and compare two recently proposed algorithms to determine the matching scales, defined in ref. [2] and [1] and referred to as BV and HMW, respectively, in what follows. In both approaches, the matching scale µ i (i = t, b, int) is determined separately for the component of the cross section involving only the top-or the bottomquark loop (µ t , µ b ), and for the top-bottom interference contribution (µ int ). The resummed results for each of these terms are then added in order to yield the best prediction for the It is worth mentioning that the integral of this equation over p ⊥ reduces it to an identity, i.e. all matching scales drop out from the equation and the correct normalization to the cross section with both the top and the bottom loop is maintained. The interference term, at variance with the first two, is not positive definite; in particular, it may vanish for a specific value of the Higgs mass. This will become relevant for the discussion in section 3.3 and 3.4. Note that due to the fact that the scales are determined separately for each component, it is possible to use them in any model with arbitrary relative strength of the couplings of the Higgs boson to the top and bottom quarks. On the other hand, the presence of any other colored particle running inside the loop would require a separate consideration. This case could be treated in the very same fashion using the methods described in this paper though.

Matching scale determinationà la HMW
The idea behind the HMW approach is the fact that (a) for p ⊥ m φ , the p ⊥ -spectrum should be well described by fixed-order perturbation theory, and (b) one would like to have an all-order result for an as large range of p ⊥ as possible. Let us first discuss condition (b). In all the approaches described above, the matching scale (Q res for AR, Q sh for MC@NLO, and h for POWHEG) can be seen as a measure up to which value of p ⊥ the all-order resummation is effective in the matched p ⊥ -distribution. Formally, one would thus like to choose the matching scale as large as possible. However, the resummation is strictly valid only in the limit p ⊥ → 0, so one cannot expect to obtain a sensible result once the matching scale gets too large. HMW therefore uses condition (a) to determine a maximum value for the matching scale. The basic idea of the HMW prescription applies in principle to any resummation/matching approach. Nevertheless, let us focus on its application to AR in what follows.

JHEP01(2016)090
Even though the matched expression of eq. (2.2) will eventually converge to the fixedorder result for p ⊥ → ∞, 17 this transition may happen only at very large p ⊥ . Typically, for large values of the matching scale, the integral over p ⊥ from zero to m φ becomes rather large, and may even overshoot the total NLO cross section. Due to the unitarity constraint of eq. (2.4), it follows that dσ/dp ⊥ can deviate significantly from the fixed-order result for p ⊥ ∼ m φ and may even turn negative in order to compensate for the excess at small p ⊥ . This spoils the whole idea behind matching the resummed with the fixed-order result and thus defines an upper limit on the matching scale.
There is certainly quite an amount of arbitrariness in this procedure, in particular: in what range should the matched result agree with the fixed-order result, and to what degree? This arbitrariness has to be taken into account in the estimate of the theoretical uncertainty. In this paper, we define the so-called HMW approach by following ref. [1], where a particular set of these parameters was defined, from which the maximal matching scales were determined.
To be precise, Q max res is defined as the maximum value of Q res for which the resummed The default matching scale Q is then defined to be half of that maximum value. As it turns out, the choice of the central matching scale as defined above indeed leads to a behavior of the matched result in the large p ⊥ region which is very close to the fixed-order result.
As pointed out above, this procedure is applied separately to the top-and the bottomquark induced contribution to the cross section, and to the top-bottom interference term, resulting in three different HMW scales Q t , Q b , and Q int , respectively.

Matching scale determinationà la BV
The validity of the resummation formalism relies on the soft and collinear factorization of the squared matrix elements describing real parton emissions. The factorization in the soft limit can be demonstrated in a straightforward way thanks to the fact that, for increasing radiation wavelengths, the details of the hard scattering process are not resolved, independently of all the other kinematic details of the emitted parton. The discussion of the collinear factorization is more complex. In the BV approach, the accuracy of the collinear approximation in the gluon-fusion process is discussed, at partonic level, in the presence of an exact description of the top and bottom quarks running in the virtual loop. The procedure to determine the scales is described by the following steps. The exact squared matrix elements of the subprocesses gg → gH and qg → qH are compared with their collinear approximation. A deviation by more than 10% from the exact result signals that the collinear approximation breaks down. 18 The upper limit w of the range of Higgs transverse momenta where the collinear approximation is accurate is chosen as the value for the matching scale in any hadron level calculation, either with analytic resummation or in a PS Monte Carlo. The two partonic subprocesses initiated by gg and by qg have a different collinear behavior, which leads to two different scales w gg and w qg ; the final scale w is computed as the JHEP01 (2016)090 average of the two previous values, weighted differentially by their relative importance to the transverse momentum distribution of the Higgs, in the p ⊥ range between w gg and w qg .
The three BV scales associated with the top, the bottom and the interference term will be denoted by w t , w b and w int , respectively.

Qualitative comparison of the two approaches
Since the matching scale (or specifically Q res , Q sh , h) is unphysical, its choice is formally arbitrary, and any prescription for its determination is necessarily heuristic. The BV and the HMW approach are complementary in at least two aspects. While BV works at the partonic level and considers the low-p ⊥ region, the HMW approach uses the large-p ⊥ region of the hadronic distribution in order to choose a value for the matching scale. Furthermore, BV does not make any reference to the specific method of resumming the logarithmic terms at small p ⊥ . The resulting scales may be interpreted as POWHEG's h, MC@NLO's Q sh , the resummation scale Q res of AR, or any other matching scale characteristic for the separation of the soft-collinear from the hard region. The scales determined through the HMW approach, on the other hand, do in principle depend on the underlying resummation technique. We will discuss this issue in the light of the three resummation codes considered here in the next section.
On the other hand, since both approaches separately treat the top-, bottom-, and the interference term, the resulting scales obtained here and in refs. [1,2] are independent of the respective Yukawa-couplings and can be applied to the SM, the 2HDM, and other models where the gluon-Higgs coupling is predominantly mediated by the third generation of quarks. In fact, the scales only depend on the CP parity and on the mass of the Higgs boson.
Considering the differences between the two approaches, it is not surprising that the numerical values of the resulting scales are different. Since the constraints that are adopted by the two groups act on different parts of the p ⊥ spectrum (low p ⊥ in the BV case, large p ⊥ in the HMW case), the spread of the results is likely to cover in a quite conservative way the ambiguities of this scale determination. The hierarchy of the Higgs and of the quark masses determines a moderate (top) or a very good (bottom) agreement between the two groups. There is one exception where the scales of BV and HMW may differ by many orders of magnitude, namely when the LO term is much smaller than the NLO term. This only happens for the interference contribution which is not required to be positive definite. Since the resummed contribution is always proportional to the LO term (apart from corrections due to the virtual contributions which are small compared to the total cross section), it will also be small in these cases, and the distribution will be given almost completely by the hard emission from the NLO term. It follows that the BV-scale will vanish as the LO term tends to zero, because the collinear approximation fails for any value of p ⊥ > 0. On the contrary, since the resummed curve is almost identical to the fixed-order one in this case, and since the HMW algorithm looks for the largest scale that fulfills the HMW criteria, the resulting matching scale will tend to be very large.

Quantitative comparison of the two approaches
After clarifying the conceptual differences between the two approaches, we can now study the actual numerical values of the matching scales. The upper plot in figure 3 shows the three scales for the t, b, and interference contributions for scalar Higgs production in the two approaches, with m φ ≤ 600 GeV in the HMW case and m φ ≤ 700 GeV in the BV case. As outlined above, while the BV scales are independent of the resummation procedure, this is not the case for HMW. The numbers shown in figure 3 are based on AR. In principle, separate HMW matching scales should be determined for MC@NLO and POWHEG. For MC@NLO, it turns out though that the HMW scales of figure 3 lead to numerical results close to the ones obtained in the AR approach, as the matching procedure is indeed rather similar in the two cases. On the contrary, POWHEG consistently exceeds the fixed-order α 3 s distribution (referred to as fNLO [56] in the following) in the high-p ⊥ tail even for very small matching scales h. Any HMW criterion which requires a transition to fNLO at large p ⊥ becomes questionable in this case. We therefore do not attempt to determine separate HMW scales for POWHEG, but simply use the ones of figure 3. Let us add that the situation changes significantly when applying the so-called mPOWHEG modification, to be defined in section 4.2.3. Similar to MC@NLO, the HMW scales of figure 3 lead to reasonable results in this case.
Concerning the top-induced contribution, w t exhibits a non-trivial structure in a broad m φ range around the top-quark threshold, between about 220 and 380 GeV. The corresponding HMW scale Q t , on the other hand, increases quite steadily with the Higgs mass, and the top-quark threshold only has a very mild effect, if any. In addition, the overall growth of w t with m φ is stronger than for Q t . Nevertheless, over the whole m φ ≤ 600 GeV region considered here, the two scales never differ by more than a factor of two. Since the resummation uncertainty in the two approaches will be estimated by varying the matching scales by a factor of two, we can expect consistent results in cases where the top contribution is dominant.
On the other hand, the BV and HMW scales for the bottom-induced contribution are in much better agreement. The BV scale w b exhibits a slightly steeper dependence on m φ than Q b , but the difference between the two remains quite small for all m φ 600 GeV. Except for very small Higgs masses, w b and Q b are considerably larger than m b . In the BV approach one observes that this result is driven mainly by the behavior of the partonic gg → φg channel [2] (see also ref. [75]), while the qg → φq channel would suggest a value for the bottom matching closer to m b [49].
For the scales of the interference contribution, both approaches lead to a very similar slope in their respective matching scales as m φ increases from about 40 to 320 GeV, even though their absolute values differ significantly. While w int remains below about 25 GeV, Q int is always larger than 20 GeV as long as m φ 30 GeV. For larger values of m φ , the two interference scales show a very different behavior: the BV scale w int slowly decreases until it reaches the value w int = 0 at about 590 GeV, which happens precisely where the interference term of the total cross section at LO vanishes. In contrast to that, the slope of the HMW scale increases beyond m φ ≈ 340 GeV, and Q int assumes its maximal value of about 80 GeV at m φ ≈ 570 GeV. A similar feature is observed around m φ = 30 GeV. The reasons for the  20  40  60  80  100  120  140  160  180  200  220  240  260  280  300  320  340  360  380  400  420  440  460  480  500  520  540  560  580  600  620  640  660  680  700 Q,w (GeV) qualitative differences in this case have already been discussed in section 3.3; they are the clearest manifestation of the different ideas behind the BV and the HMW method.
The corresponding plot for a CP-odd Higgs boson is shown in figure 3; the general behavior for all the three contributions is quite similar to the scalar case and requires no separate discussion. In this section, we present a quantitative comparison of the impact of the BV and of the HMW scale determinations, using three different codes: MoRe-SusHi, which implements the analytic resummation of the p ⊥ spectrum, gg H 2HDM and aMCSusHi, which apply the POWHEG and the MC@NLO method, respectively. Although the theoretical basis underlying these three codes is the same, namely resummation based on soft and collinear factorization, the specific algorithms involved are quite different, see section 2. Therefore, we will try to disentangle effects due to these different implementations from those arising from the different matching scales.
The uncertainty band due to a variation of the matching scale is obtained by the following procedure: given a set of reference values (µ t , µ b , µ int ) for the three matching scales according to eq. (3.1) (adopting either the BV or the HMW approach) we consider all the possible combinations which can be generated by taking half and twice these values, or the reference values themselves; for each setting, we compute the transverse momentum distribution; collecting all these results, we take the envelope for each p ⊥ value, i.e. the minimum and the maximum values among all the simulations. 19 For the analytic resummation, which consistently matches the fixed-order results at large transverse momenta, we follow ref. [1] and apply an additional factor to the error band which damps it towards large values of p ⊥ . This takes into account condition (a) from section 3.1, according to which resummation should not have a big impact on the large-p ⊥ region, and thus also not on its theoretical ambiguities. As we will see, the shape of the uncertainty band as derived from the variation of the matching scales has a feature common to all the codes. In a region just above the peak of the p ⊥ distribution, the band is relatively narrow. This is a consequence of the unitarity constraint, which establishes an anticorrelation between the low-and the high-p ⊥ tails of the distribution. The precise position of this region depends on the position of the peak of the resummed distribution, on the total variation of the cross section, as well as on the central value of the matching scale; the interplay of these factors determines the precise shape of the uncertainty band.

Setup and representative scenarios in the 2HDM
This section defines the phenomenological scenarios considered in this paper. They have been designed to highlight possible interplays between the top-quark and the bottom-quark mediated amplitudes.
We start with the SM, where the p ⊥ spectrum of the Higgs boson is known through NLO+NLL including the full dependence on the quark masses [48] and through NNLO+NNLL for the top-quark induced terms and by assuming the limit m t → ∞ [40,58,19 Recall that the accuracy of AR is NLL, while the NLO+PS MCs only resum consistently all LL terms and partially the NLL ones. Therefore, the higher order terms probed by the scale variation procedure are slightly different in the two cases.  59]. 20 In the present study, for uniformity with the BSM codes where only NLO QCD accuracy on the total cross section and NLO+(N)LL results on the Higgs transverse momentum distribution are available, we restrict also the SM analysis to this level of accuracy.
While there is hardly any controversy that the "characteristic scale" for the top-quark induced contributions in the SM should be of the order of the Higgs boson mass (m h = 125 GeV), already in this case the BV and the HMW methods provides us with a more quantitative result for the matching scale which turns out to be close to the often adopted choice of m h /2 in the SM, but becomes significantly smaller than m φ towards larger values of m φ (see figure 3). One of the main subjects of this paper, however, is the question of how to take into account the bottom-quark induced contribution. In the SM, the effect of the bottom quark is suppressed by the Yukawa coupling, and therefore differences in this treatment have limited effect on the overall momentum distribution. In models with an extended Higgs sector, however, the bottom-quark Yukawa coupling can be significantly enhanced, at least for some of the Higgs bosons of such theories.
One of the simplest extensions of the SM in this respect is the 2HDM. Therefore, we focus on 2HDM scenarios which we devise in order to enhance specific contributions to the cross section. 21 The conclusions of our study, however, trivially generalize to other models where the gluon-Higgs coupling is predominantly mediated by top and bottom quarks (e. g. most of the experimentally viable parameter space of the MSSM). Since the aim of this paper is a conceptual one, we disregard any phenomenological constraints on the 2HDM parameter space, in general; we do not consider the case of a light Higgs boson of mass m h = 125 GeV with enhanced bottom-quark Yukawa coupling though, due to the obvious conflict with experimental observations. We do, however, respect the theoretical constraints due to unitarity and triviality of the theory, as well as stability of the physical vacuum. We check these constraints with the help of the program 2HDMC [76,77]. In all scenarios, except for the low-m A scenario to be introduced at the end of this section, we set the mass of the two CP-even Higgs bosons to m h = 125 GeV and m H = 300 GeV, respectively, while the mass of the CP-odd Higgs boson is set to m A = 270 GeV.
The first scenario we consider is "scenario B" as defined in ref. [78]. For consistency with the notation in the rest of this paper, however, we will refer to it as "large-b scenario" in what follows. Since it induces a SM-like light Higgs boson, we only study the production of the heavy and the CP-odd Higgs boson for that scenario, both of which have strongly enhanced couplings to the bottom quark.
As a modification of this, we use the same parameters as in the large-b scenario, except that we set tan β = 1. This will be referred to as "large-t scenario" in what follows. It is designed to result in a top-quark dominated cross section for the heavy and pseudo-scalar Higgs boson. Again, the light Higgs is very SM-like and will not be considered in this scenario. 20 Recall that the (N)NLO accuracy denoted here refers to the underlying total cross section (including terms up to α 3 s and α 4 s , respectively) and corresponds to (N)LO accurate predictions at large transverse momenta. 21 More precisely, we only consider the type-II 2HDM, in which one doublet generates the masses of the up-type quarks and the other of down-type quarks and charged leptons. In any case our results are directly applicable also to the other 2HDM types. Finally, we devise a set of rather pathological scenarios where the LO cross section receives a large top-bottom interference contribution. The parameters of these "large-int scenarios" have to be chosen differently for each of the neutral Higgs bosons.
The precise definition of all scenarios, together with the top-, bottom-and interference component of the total inclusive cross sections at LO and NLO, is given in table 1. Note that, while for the light and the pseudo-scalar Higgs the absolute value of the interference term in the large-int scenarios amounts to more than 100% of the total cross section, we did not manage to find a parameter point for the heavy Higgs which has a similarly large interference term while still respecting the theoretical constraints of unitarity, stability, and perturbativity.
Let us emphasize again that most of these scenarios are in vast conflict with experimental observations; they only serve as theoretical benchmarks for the study of resummation ambiguities in the Higgs p ⊥ distribution. For phenomenologically viable 2HDM benchmark points we refer the reader to ref. [79].
We further investigate one phenomenologically interesting scenario with a very low pseudo-scalar Higgs boson mass. In this case we chose a scenario of ref. [80] that meets all theoretical as well as experimental constraints. This scenario is referred to as low-m A in table 1. The masses of the three Higgs bosons are m h = 125.5 GeV, m H = 507 GeV and m A = 29.9 GeV. For such a low Higgs-boson mass, the gluon-fusion process is particularly important, since its cross section is highly enhanced with respect to Higgs production in association with bottom quarks 22 and dominant even at large tan β [80].
All the numerical results are computed for the LHC, with a center-of-mass energy of √ S = 13 TeV. We use the MSTW2008nlo68cl PDF set [83] through the LHAPDF6 library [84]

Numerical results in the different scenarios
This section contains the results obtained with the three codes under study, in the various scenarios defined in section 4.1. Since our focus is on the uncertainties inherent to the matching procedure, we compute the uncertainty band by varying only the matching scales, as described in section 4.
For each scenario we performed two different studies: 1. The predictions for AR, MC@NLO and POWHEG are compared by using the same numerical values for the respective matching scales; in this way we can assess the impact of the higher-order QCD terms which are included in different ways in the three codes, due to the different matching procedures.
2. For each code, the predictions obtained by setting the matching scale to the BV and to the HMW values are compared; this allows us to assess the sensitivity of each code to a matching scale variation, the impact of the BV vs. the HMW prescription, and the overlap of the respective uncertainty bands.

Results in the Standard Model
Even though for the SM one may expect consistency among the three codes, it will be instructive to start our discussion with this case, because it highlights certain generic features of the individual approaches which will be carried forward also to some of the other scenarios considered in this paper. Figure 4 shows the shape of the transverse-momentum distribution (i.e., the integral of each curve is normalized to one) for a SM Higgs of m h = 125 GeV in the range 0 ≤ p ⊥ /GeV ≤ 400. In the two upper plots, we compare the results of the three codes, setting the matching scales to the same numerical values: BV scales in the left and HMW scales in the right panel (cf. item i) above). Each of the lower plots, on the other hand, was obtained with one particular code (left: AR; center: MC@NLO; right: POWHEG); the different lines correspond to different values of the matching scales, BV and HMW (cf. item ii) above). All plots contain the fNLO result as a reference. In order to facilitate the discussion, we show the same plots in figure 5, but with enlarged low-p ⊥ region (0 ≤ p ⊥ /GeV ≤ 100). Apart from the fact that, as expected, the three codes yield compatible results within uncertainties (at least for p ⊥ ≤ m h ), it is worth noting the following characteristics: • The central POWHEG and the MC@NLO spectra are in excellent agreement (at the fewpercent level) between about 10 < p ⊥ /GeV < 130, while they differ by about 20% from the central AR prediction in most of this region.
• The peak position for both POWHEG and MC@NLO is at around p ⊥ = 12 GeV, while the peak for AR is about 2 GeV below that.

JHEP01(2016)090
• The central value of AR approaches fNLO at the level of 5% above p ⊥ ≈ 130 GeV; for MC@NLO, such a transition to the fixed-order curve occurs at about p ⊥ ≈ 180 GeV, while POWHEG always remains about 20% above fNLO in the tail of the distribution. The latter characteristic is a general feature of the default POWHEG matching, as described in section 2.2.2, and will be analyzed in detail in section 4.2.3.
• Above p ⊥ ≈ 130 GeV, the uncertainty band for AR is suppressed due to the damping factor introduced in eq. (4.1). Concerning the two MCs, in MC@NLO the uncertainty, expressed in unit of the AR central value, is of the order of ±10% from p ⊥ = 130 GeV all the way up to p ⊥ = 400 GeV, while for POWHEG it decreases uniformly from ±20% to ±10%.
• The MC@NLO and POWHEG bands nicely overlap for p ⊥ ≤ 100 GeV; above this value the overlap is only partial, because of the different central predictions.
• Towards smaller values of p ⊥ , the uncertainty bands for all three codes develop a bulgy structure with a maximum of about ±20% (±35%) for AR (the Monte Carlos) and a minimum of a few percent slightly above the peak position.
• Towards even smaller values of p ⊥ , the AR uncertainty band quickly grows to the 100% level, and the POWHEG band to about ±40%. Only the MC@NLO band increases to a moderate ±15%.
Since, for a SM Higgs of m h = 125 GeV, the BV and HMW scales are quite close to each other (see figure 3), the two upper plots in figure 4 are very similar. In order to see this more explicitly, we study the impact of the different scale choices (BV or HMW) within one particular code (left: AR; right: MC@NLO) in the lower three plots of figure 4. From the latter we can indeed see that the results obtained in the two cases are in very good agreement with each other, both with respect to the central value and the uncertainty band; the only difference being at very low transverse momenta (p ⊥ < 10 GeV) in the central predictions of AR and POWHEG. This can be traced back to the different matching scales for the interference term, which constitutes the largest contribution induced by the bottom-quark loop in the SM.

Results in the large-t scenario
Let us now consider the production of a heavy Higgs with m H = 300 GeV in the large-t scenario, where, similar to the SM, the bottom-loop and the interference contribution play only a minor role. In contrast to the SM, however, the matching scales for the top-loop contribution derived using BV and HMW differ by almost a factor of two. The results for the p ⊥ distribution are shown in figure 6.
Using BV scales, one notices that AR deviates quite significantly (∼ 50%) already at p ⊥ = 400 GeV from the fNLO result; this deviation tends to further increase towards larger p ⊥ values. This is not unexpected since the HMW scales are designed to guarantee similarity between the resummed and the fNLO curve at large p ⊥ . Scale choices larger than the values determined by HMW will therefore necessarily lead to a deviation from the fNLO in that  Figure 5. Same as figure 4, but restricted to the low-p ⊥ region.
region. In the case of POWHEG and MC@NLO, their canonical high-p ⊥ behavior starts to appear only for relatively large values of p ⊥ , the reason for this being again the relatively high value of the scale for the top contribution (w t = 111 GeV). Indeed, the agreement between the two Monte Carlos turns out to be excellent, at least up to p ⊥ values as large as the Higgs mass. Despite the large deviations of AR in the tail and the much softer AR spectrum, all approaches are compatible within uncertainties at small to intermediate transverse momenta (p ⊥ 200 GeV). It should be noted that this is partly due to the fact that the uncertainty bands are significantly larger (almost by a factor of two) than in the SM. Using HMW scales, the transition to the high-p ⊥ region is more similar to the SM case (including the consistent overshooting of the POWHEG spectrum) because of the relatively flat dependence on the Higgs mass of the HMW scales, implying Q large-t t /Q SM t 1.3. The results from the three codes appear to be more compatible, in particular in the tail of the distributions, mainly due to the different AR behavior. The bulges of the uncertainty bands above the peak position extend to considerably larger values of p ⊥ , and their width is similar to what we observed for the SM.   Despite the apparent differences between the left and the right upper plots, the dedicated analysis of the impact of the scales choice in the lower plots reveals that the results are nicely compatible within the respective uncertainty bands. The only exception from this occurs for AR at large p ⊥ . But one should bear in mind that the uncertainty band for AR is manually suppressed at large p ⊥ , a procedure that should strictly be applied only when HMW scales are used. The observation that the width of the BV-bands is larger than for HMW is due to the fact that w t is more than twice as large than Q t , and thus the variation by a factor of two has a bigger impact on the final result.
The results for the CP-odd Higgs boson in the large-t scenario are very similar to the ones of the heavy CP-even Higgs shown here. We therefore refrain from showing them here.

Results in the large-b scenario
The large-b scenario produces large bottom-Yukawa couplings for the heavy and the pseudoscalar Higgs boson which renders the bottom-loop induced contribution by far dominant. Since the associated matching scales w b and Q b are very close to each other, any difference in the p ⊥ distributions are due to the conceptional variants of the matching in the three codes under consideration. Figure 7 shows the results for a heavy Higgs boson using HMW scales (the corresponding plots for BV scales are identical for all practical purposes).
Let us first discuss the upper left plot of figure 7, where the curves are displayed in the same way as in the corresponding plots of figure 4 and figure 6 for the SM and the large-t scenario, respectively. While the large-p ⊥ behavior for AR and MC@NLO is similar to the large-t scenario, POWHEG produces a spectrum that is significantly harder, exceeding the fNLO result by about 50% for p ⊥ 200 GeV. Apparently, the specific matching procedure of POWHEG has a significant impact on the large p ⊥ region, where the parton shower, based on the soft/collinear approximation, is outside its region of validity. Between 10 p ⊥ /GeV 130, on the other hand, POWHEG and AR agree within 10%, while MC@NLO is significantly higher in that region. The central predictions of the two Monte Carlo results are in reasonable agreement only at small transverse momenta (p ⊥ 30 GeV). Moreover, the size of the error bands is very different in the two Monte Carlo approaches: the MC@NLO band blows up to O(100%) around p ⊥ ∼ 125 GeV; the POWHEG band remains very narrow over the whole range, indicating an uncertainty even smaller than in the SM. The fact that all approaches lead to compatible predictions below p ⊥ ≈ 200 GeV is mainly due to the large MC@NLO band.
Since the large-b scenario reveals the differences between the three codes under study in the most striking way, it will be instructive to investigate them in more detail within this scenario. Consider first the POWHEG approach. As it turns out, both observed features -enhanced high-p ⊥ tail and small uncertainty band -can be tackled by the same modification of the matching procedure: in the original POWHEG approach, the scale t 1 (see section 2.2) for each event is identified with the transverse momentum of the first emission. If the latter is very large, the shower will act up to scales which are way beyond the validity range of the underlying approximations. If instead we restrict t 1 for all remnant events (the R f -term in eq. (2.6)) to remain below the matching scale (e.g. BV or HMW), one obtains the result shown in the lower left plot of figure 7 (magenta, solid curve with stars). Since the high-p ⊥ tail of the distribution is driven by the remnant events, the above restriction of t 1 ensures a transition to the fNLO curve. We will refer to this modified approach as "mPOWHEG" in what follows. 23 It happens that also the uncertainty band is more similar to the other approaches for the mPOWHEG result (see again the lower left plot of figure 7).
The fact that the mPOWHEG modification applies only to the remnant events ensures that the formal accuracy of the original POWHEG approach remains unaffected. In this respect, it is reassuring that at low transverse momenta (p ⊥ 100 GeV), the modifications have practically no effect (compare the POWHEG with the mPOWHEG curve in figure 7), as this region is controlled by the soft events (the first term in eq. (2.6)), which remain unchanged. 23 Modifications of the shower scale setting were studied also in ref. [86] in order to improve the theoretical prediction for dijet production.

JHEP01(2016)090
We checked that the fNLO-like high-p ⊥ behavior of mPOWHEG is a generic feature and not specific to the underlying phenomenological scenario or the details of the parton shower. A more quantitative study of its numerical impact has to be deferred to a future study though. Stressing again that, to our understanding, mPOWHEG is only a logarithmically subleading modification of POWHEG (or rather its interface with the parton shower), one may consider it as a viable alternative whenever an fNLO-like high-p ⊥ behavior is required.
Let us now discuss the MC@NLO approach in this scenario, featuring a peculiarly large uncertainty regarding shower scale variations, which in turn leads to very different shapes of both the central predictions and the error band with respect to the other two approaches. Recall that in the MC@NLO implementation of MadGraph5 aMC@NLO, the values for the shower scale of the soft events follow a specific distribution with a peak at the matching scale µ i (i = t, b, int; e.g. BV or HMW), see section 2.2.1. We find that, in the large-b scenario, a restriction of the range of that distribution has a significant effect on the central MC@NLO prediction, which is not surprising given the large associated uncertainty. In the limit where this distribution turns into a δ-function δ(Q sh − µ i ), also the size of the uncertainty band is strongly reduced, as can be seen in the lower right plot of figure 7 (orange, solid curve with full boxes). We checked that, besides a slightly earlier matching to the fNLO prediction controlled by the remnant events, the observed features are due to the restricted shower scale range accessible to the soft events. This study simply shows that the predictions in bottom-quark dominated scenarios depend strongly on details of the matching procedure and thus should be attributed with large uncertainties.
Finally, the upper right plot of figure 7 shows the comparison of the results obtained with the modifications done for MC@NLO and POWHEG. Indeed, the modified predictions turn out to be in much better agreement, in terms of both the central curves and the uncertainty bands.
Similar to the large-t scenario, the distributions for the pseudo-scalar Higgs in the large-b scenario largely resemble the ones of the heavy Higgs shown in figure 7, and we do not need to discuss them separately at this point.

Results in the large-int scenario
The matching scales of the interference term in the BV and the HMW approach exhibit quite a different behavior, as shown in figure 3. It will thus be interesting to see the distributions in the various approaches in a scenario with a particularly large interference term. Note, however, that the interference term in the large-int scenarios defined in table 1 is always negative and competes with a similarly large top and/or bottom contribution. Our 2HDM parameter scan did not reveal a truly interference-dominated scenario, where both the top and the bottom contribution are small compared to the interference term, and which passes the unitarity, stability, and perturbativity checks. Figure 8 shows the results for the light Higgs boson, where the modulus of the top, the bottom, the interference term, and the total cross section each are roughly of the same size (∼ 2 pb, see  the BV and the HMW approach for the light Higgs, while the interference matching scale is about a factor of three larger in the HMW approach. Using HMW scales, the comparison of the distributions in the various approaches leads to a picture that is quite similar to the one that we found in the SM: good agreement among the Monte Carlos, softer spectrum of AR, compatible results within uncertainties of all approaches (at least for p ⊥ ≤ m h ). The error bands, on the other hand, are significantly larger than in the SM. They amount up to about 24 ±60% for AR and MC@NLO, and about ±40% for POWHEG.
Using BV scales, on the other hand, the width of all uncertainty bands is strongly reduced (up to about ±20% and ±30% for AR and the Monte Carlos, respectively), the 24 We disregard the region very close to the threshold (p ⊥ = 0), where, as stated before, the AR uncertainty band blows up.   reason being again the fact that a large Q int (as in HMW) also induces a large interval [Q int /2, 2Q int ]. While also for BV the results at small p ⊥ are compatible within uncertainties and the high-p ⊥ behavior is very reminiscent of the one observed in the SM, there are some substantial differences at very small transverse momenta (p ⊥ < 40 GeV) between the shapes predicted by the two Monte Carlos. For better visibility, figure 9 shows in the left panel the BV plot with an enlarged low-p ⊥ region. Clearly, the MC@NLO spectrum is harder, while the POWHEG one actually gets negative in the first bin. 25 Since the interference term is not positive definite, due to its definition by subtraction, the distribution may turn negative in scenarios where its contribution is large. Clearly, this behavior appears to be amplified at small p ⊥ if the matching scale is particularly small. At least this specific case, therefore, leads to a similar conclusion as pointed out in ref. [54], that very low scales are not well suited for NLO matched parton shower predictions. Indeed, the behavior just observed is even enhanced in the large-int scenario of the pseudo-scalar Higgs with m A = 270 GeV, see figure 9 (right panel). The MC@NLO distribution displays some odd behavior between about 10 < p ⊥ /GeV < 60, developing an almost linear behavior rather than the ordinary Sudakov shoulder. The POWHEG prediction, on the other hand, becomes negative in the first and the second bin of the distribution (i.e., below p ⊥ = 10 GeV). There are mainly two reasons for these features being intensified in the large-int scenario of the pseudo-scalar Higgs: first, the ratio of the interference matching scale to the Higgs mass is even smaller than for the light Higgs; second, the relative contribution of the interference term is larger than in the light Higgs case. The latter can be inferred from table 1: for the pseudo-scalar Higgs, the large-int scenario corresponds to 25 Note also that the AR curve for BV scales turns negative at very low p ⊥ .   an interference term whose modulus is about 30% larger than the top, 50% larger than the bottom contribution, and 270% larger than the total cross section. The other results in the large-int scenario of the pseudo-scalar Higgs are shown in figure 10. The conclusions are essentially the same as the ones for the light Higgs, which will thus not be repeated here. It is worth mentioning though that the relative size of the uncertainty bands for both BV and HMW scales are increased with respect to the light Higgs case, the reason being again the larger value of the matching scale and the thereby enlarged variation range.

Results in the low-m A scenario
Also the low-m A scenario produces large bottom-Yukawa couplings, but in this case the Higgs and the bottom mass are significantly closer than in the large-b scenario. Thus, all contributions besides the bottom-loop induced one are negligible and Q b is quite close  to w b . Figure 11 shows, therefore, only the comparison among the three codes for the BV scales, while the corresponding HMW results are identical for all practical purposes.
Looking first at the left plot in figure 11, we notice that the characteristic features observed in the large-b scenario appear to be amplified in this case: the MC@NLO uncertainty band blows up to O(150%) around p ⊥ ∼ 25 GeV; the POWHEG curve exceeds the fNLO result by about 100% for p ⊥ 50 GeV with a very small uncertainty band. The central curves are vastly different between the two Monte Carlo approaches except for the first two bins.
In the right plot of figure 11 we show the corresponding curves for the modified approaches (MC@NLO with fixed shower scale, mPOWHEG). Evidently, the agreement among the results is significantly improved. Also in this case, a restriction of the shower scale in POWHEG closes the large gap to fNLO in the tail of the distribution, such that the mPOWHEG curve nicely approaches fNLO at large p ⊥ . The MC@NLO uncertainty band is very similar to the one of AR and the central predictions of the Monte Carlos are almost on top of each other. Note that for such small Higgs boson masses, the shower scale distribution as defined in section 2.2.1 ranges down to very low values (Q sh ∼ 1 GeV) which are not well suited for a Monte Carlo approach, both from a physical and a technical point of view. This can be avoided by using narrower distributions -a δ-function in our case -for Q sh , thus leading to a better agreement with the other codes and a more reasonable size of the uncertainty band. 26 26 In the default implementation of MadGraph5 aMC@NLO [56], there is a hard cut Q sh ≥ 3 GeV that avoids such low values; for this study, however, we have removed this cut.

Conclusions
A reliable theoretical prediction of the Higgs transverse-momentum distribution must include an estimate of the uncertainties associated with the matching of fixed-and all-order results. The latter are relevant to describe the large-and the low-p ⊥ parts of the spectrum, respectively. The matching procedure is not unique; in this paper, we studied three different methodological approaches and their dependence on the respective matching scale.
In general, the gluon-fusion process is characterized by two external, physical scales: the mass of the Higgs boson and the mass of the quark running in the loop of the gluongluon-Higgs vertex. Since in the SM, where the top quark is dominant, these two scales are quite close to each other, the issue of setting the matching scale is less problematic. In BSM scenarios with enhanced bottom Yukawa coupling, however, the proper choice of the matching scale has become a matter of debate.
In this paper we presented two distinct comparisons: 1. The determination of a suitable central value for the matching scales associated to the top, bottom, and interference terms, as proposed in the HMW [1] and BV [2] approaches, is compared qualitatively and quantitatively. Figure 3 summarizes the main outcome as a function of the Higgs mass, both for a CP-even and a CP-odd Higgs boson.
2. The predictions of three resummation codes to the problem of the matching (AR, MC@NLO and POWHEG), all at NLO QCD accuracy concerning their normalization, but different in their logarithmic accuracy, are compared at the level of their central values and of the respective uncertainty bands, obtained through variation of the matching parameter by a factor of two around a central value. The latter has been fixed following the prescriptions of point i). In this comparison we considered the production of a Higgs boson in the framework of the SM and a number of 2HDM scenarios which are representative of the different possible interplays between topand bottom-quark effects in the gluon-fusion scattering amplitude.
In these comparisons, it is important to keep in mind that resummation considerably improves the theoretical prediction in the low-p ⊥ region, while one cannot expect the resummed (and matched) result to be in any way superior to the fixed-order prediction at intermediate and high-p ⊥ . We showed that the discrepancies among the codes outside the low-p ⊥ region can be greatly reduced by a modification of the treatment of formally subleading logarithmic terms.
In the low-p ⊥ region, we find very compatible results among the three codes in all scenarios, independent of whether one uses BV or HMW scales. Even in cases where the central scales suggested by the two methods differ relatively strongly, the respective p ⊥ distributions are consistent with each other once their error band is taken into account. For the interference term, this reflects the complementary reasoning behind the BV and HMW methods which allows them to diverge as long as the sensitivity of the p ⊥ distribution on the matching scale is small.

JHEP01(2016)090
Let us summarize the outcome of our study as follows: • The matching uncertainty on the prediction of the Higgs transverse momentum distribution amounts to one up to several tens of percent, depending on the specific value of p ⊥ and on the model under study. This uncertainty should be combined with the usual perturbative uncertainty estimated via renormalization and factorization scale variations.
• As far as small transverse momenta are concerned, we find reasonable agreement among the codes irrespective of the specific matching scale choice (BV or HMW), although the central AR result shows a generally much softer spectrum as compared to the central prediction of the Monte Carlo codes. The agreement among the three codes is in part caused by an increase of their error band towards very low p ⊥ (essentially in the first 5 GeV bin of the distribution), which is most apparent though for the AR result. Note, however, that the cross section is typically very small in this region. Let us recall that these results have been obtained with the Pythia8 parton shower; variations with the parton shower are left for future studies.
• In the intermediate p ⊥ range, the estimate of the width of the matching uncertainty bands shows differences in the three approaches under study, which are more pronounced in the bottom dominated scenarios. In this last case the size and shape of the uncertainty bands strongly depends on the matching formulation, as it has been explicitly shown in the MC@NLO case with the use of different probability functions to extract the value of the shower scale.
• In the large p ⊥ range, where all the codes have only LO accuracy, sizable differences appear in their predictions, where the POWHEG result is systematically larger than the MC@NLO and AR curves, which in turn are compatible with the fNLO one. One source of this discrepancy has been identified in the different treatment of high-p ⊥ multiple parton emissions beyond the first one, still allowed in POWHEG and suppressed in the other two cases. We expect the precise size of the discrepancy to depend on the details of the parton shower; quantitative studies will be deferred to a future publication. The codes compared in this note provide either a transition to the LO prediction at large p ⊥ (AR and MC@NLO), which has a well defined perturbative accuracy, but misses the effects of additional radiation, or a prediction that includes multiple parton emissions at all transverse momenta, but describes them by means of a PS (POWHEG), which is not adequate at large p ⊥ because it is based on the soft/collinear approximation. It was shown that a simple modification of the interface between POWHEG and the PS (referred to as mPOWHEG in the text), which restricts the shower scale of the remnant events, can be applied to ensure a consistent matching to the fNLO prediction in the tail of the distribution, similar to the other two approaches.
The approaches studied here do not lead to identical results; not only their central values, but also their matching-induced uncertainty bands show characteristic differences.

JHEP01(2016)090
Nevertheless, it is fair to say that, in the low-p ⊥ region and including these error estimates, the results obtained by following either the BV or the HMW approach, and by employing either of the three codes are consistent with each other. This allows the conclusion that any of these approaches leads to a valid prediction for the low-p ⊥ distribution once the matching as well as the perturbative uncertainties are taken into account. The former should be estimated by a variation of the respective matching scale around its central value, as given by BV or HMW, the latter by a variation of the renormalization and factorization scales. In the intermediate-and high-p ⊥ region, where large PS effects become quite doubtful, all approaches can be made compatible with the fNLO result by suitable minor modifications.
Let us point out that the new generation of Monte Carlo event generators, merged at NLO+PS [87][88][89] or NNLO QCD accurate for quantities inclusive over the radiation [60][61][62][63][64], with NLO QCD accuracy in the description of the Higgs transverse momentum spectrum, offers a more accurate description of the p ⊥ spectrum, however only in the SM-like scenario. The evaluation of the associated matching uncertainties is left to a future study.  Table 3. Table of the matching scales (in GeV) in the HMW and BV approach for a CP-odd Higgs boson. A dash is used to indicate the case where the determination procedure of a scale has not been successful.