Modeling BSM effects on the Higgs transverse-momentum spectrum in an EFT approach

We consider the transverse-momentum distribution of a Higgs boson produced through gluon fusion in hadron collisions. At small transverse momenta, the large logarithmic terms are resummed up to next-to-leading-logarithmic (NLL) accuracy. The resummed computation is consistently matched to the next-to-leading-order (NLO) result valid at large transverse momenta. The ensuing Standard Model prediction is supplemented by possible new-physics effects parametrised through three dimension-six operators related to the modification of the top and bottom Yukawa couplings, and to the inclusion of a point-like Higgs-gluon coupling, respectively. We present resummed transverse-momentum spectra including the effect of these operators at NLL+NLO accuracy and study their impact on the shape of the distribution. We find that such modifications, while affecting the total rate within the current uncertainties, can lead to significant distortions of the spectrum. The proper parametrization of such effects becomes increasingly important for experimental analyses in Run II of the LHC.


JHEP03(2017)115
known since long time [55,56] including the full quark-mass dependence. It took nearly ten years until the O(α 4 S ) corrections were computed [57][58][59][60]. These were carried out in the heavy-top limit (HTL, i.e. m 2 t M 2 H , p 2 T H ). Finite top-mass effects on the Higgs p T distribution at O(α 4 S ) were estimated in refs. [61][62][63]. Recently, results on Higgs+jet production at O(α 5 S ) were also obtained in the HTL [64][65][66]. In the low-p T region (p T m H ), the convergence of the perturbative expansion is spoiled by the presence of large logarithmic terms of the form α n S ln m (m 2 H /p 2 T ). In order to obtain reliable predictions also in this region, the large logarithmic terms must be resummed to all orders [67][68][69][70][71]. It is then essential to consistently match the resummed and fixed-order calculations in the intermediate p T region, so as to obtain accurate predictions in the entire region of transverse momenta. In the case of the Higgs p T spectrum the resummation has been performed up to next-to-next-to-leading logarithmic accuracy (NNLL) and matched to the fixed-order NNLO result up to O(α 4 S ) in the HTL [70]. Finite quark-mass effects have been included in the resummed spectrum up to NLL+NLO [72,73]. The recent computation of the O(α 5 S ) corrections at high-p T , together with new available information on the logarithmic structure at the same order [74] would in principle allow to extend the accuracy of this calculation. 1 To consistently introduce deviations from the SM Higgs sector explicit models beyond the SM (BSM) can be directly studied. A complementary bottom-up framework is offered by the Standard Model Effective Field Theory (SMEFT). In such approach the Standard Model Lagrangian is extended by the inclusion of operators of higher dimension (in first approximation: dimension six), built from Standard Model fields and suppressed by powers of the New Physics scale (Λ) [76][77][78][79]. This enables a theoretically consistent parametrisation of model-independent effects of high-scale New Physics, which manifest themselves through small deviations from the SM predictions. Many groups translated the data already collected by the LHC as well as earlier experiments into bounds on the Wilson coefficients of dimension-six SMEFT operators (see e.g., refs. [80][81][82][83][84][85]).
On the other side, a significant effort has been devoted to supplement the tools used in the modelling of LHC data with the effects of appropriate dimension-six SMEFT operators (see e.g. [86][87][88][89][90][91]). This is highly relevant since in this way indirect BSM effects can be directly tested in the experimental analyses. The precision reached by the current experiments will call for theoretical improvements and effects from SMEFT operators beyond leading order [92]. Despite conflicting approaches followed in the literature, SMEFT effects should be evaluated by including all possible operators contributing to the observable (at a given order). Results for the total Higgs production cross section including modified top and bottom Yukawa couplings and an additional direct Hgg interaction have been obtained at NNLO in ref. [93] and at N 3 LO in refs. [94,95]. Studies of the prospects of future LHC runs for the determination of Wilson coefficients were performed with the use of such tools in refs. [96,97].
The inclusion of dimension-six and dimension-eight operators in the p T -spectrum has been considered in refs. [97][98][99][100] and [101,102], respectively. Strategies for extracting information on the Higgs-gluon couplings from the measurements were studied in ref. [97]. 1 Work in this direction has been recently presented in ref. [75].

JHEP03(2017)115
Most of the above studies, however, are limited to the high-p T region of the spectrum, and do not include small-p T resummation. In this paper we present a computation of the resummed p T -spectrum at NLL+NLO accuracy, with the inclusion of a set of dimension-six operators relevant for Higgs boson production.
The paper is organised as follows. In section 2 we discuss the effects on the Higgs production cross section from the inclusion of the dimension-six operators and we explicitly evaluate the modifications of the inclusive LO cross section. In section 3 we outline the computation of the p T spectrum of the Higgs boson, review the formalism of transversemomentum resummation required at small p T and describe our NLL+NLO calculation. In section 4 we present our results for the p T spectrum and study its sensitivity to BSM effects of the dimension-six operators. In section 5 we summarize our results and provide some concluding remarks.

Effective operators and their impact on the Higgs production cross section
We consider the effective Lagrangian where the SM is supplemented by the inclusion of a set of dimension-six operators describing new physics effects at a scale Λ well above the electroweak scale. In our study we consider the following four operators These operators, in the case of single Higgs production, may be expanded as: The operator O 1 corresponds to a contact interaction between the Higgs boson and gluons with the same structure as in the heavy-top limit of the SM. The operators O 2 and O 3 describe modifications of the top and bottom Yukawa couplings. The operator O 4 is the chromomagnetic dipole-moment operator, which modifies the interactions between the gluons and the top quark 2 (here σ µν = i 2 [γ µ , γ ν ]). In our convention, based on the SILH basis [103,104], we express the Wilson coefficients as factors in the canonically normalized Lagrangian.   The coefficients c t , c b and c g can be probed in Higgs boson processes. In particular, c t (and c b ) may be measured in the ttH (and bbH) production modes. 3 The coefficient c b can also be accessed through the decay H → bb. The coefficient c tg , instead, is constrained by top pair production [115].
We now consider the contribution of the effective operators in eqs. (2.4), (2.5) and (2.7) on the production cross section, while omitting, for simplicity, the bottom contribution in eq. (2.6). The relevant Feynman diagrams are displayed in figure 1. The corresponding amplitude can be cast into the form where τ = 4m 2 t /m 2 H and 1 and 2 are the polarization vectors of the incoming gluons. The contribution of the chromomagnetic operator to the function F (τ ) has been addressed in the literature with contradicting results [116,117] (see also ref. [118]). In ref. [116] it is found that the UV divergences in the bubble and triangle contributions cancel out. In the revised version of ref. [117] it is instead stated that the UV divergence is present, and it has to be reabsorbed into the coefficient c g .
Our results are consistent with the latter statement. We find with the functions (2.14) 3 See refs. [105][106][107][108] and refs. [109][110][111][112][113][114], respectively, and references therein. The 1/ divergence can be reabsorbed in the MS renormalization of the coefficient c g : where µ R denotes the renormalization scale of c g . The final result reads where In the HTL m 2 t m 2 H we have In the SM we have c t = 1 and c g = c tg = 0, so that F (τ ) → F 1 (τ ). In ref. [115] data on top production are used to extract constraints on c tg . The resulting region of allowed values of c tg has been found to be − 0.04 c tg 0.04 .
The impact on the total cross section is less than 20%. We conclude that, although smaller than the impact of c g , the effect of c tg can still be important. We note, however, that the chromomagnetic operator provides a contribution which is formally O(λ 2 t ) with respect to the others. In a strict expansion in α S it can be neglected. This is what we will do in the next section. Focusing on the impact of c t and c g , we note that the total cross section alone does not allow us to disentangle the coefficients c g and c t : (2.20) As already noted in the literature [99], the transverse momentum spectrum allows us to break this degeneracy.

Transverse-momentum spectrum
We consider the inclusive hard-scattering process where the colliding hadrons h 1 and h 2 with momenta p 1 and p 2 produce the Higgs boson H with transverse momentum p T accompanied by an arbitrary and undetected final state X. According to the QCD factorization theorem the transverse-momentum cross section is evaluated as F ) are the parton densities of the colliding hadrons at the factorization scale µ 2 F , dσ H,a 1 a 2 /dp 2 T is the partonic cross section,ŝ = x 1 x 2 s is the partonic centre-ofmass energy, and µ R is the renormalization scale. 4 In the low-p T region (p T m H ), the perturbative expansion is affected by large logarithmic terms of the form α n S ln m (m 2 H /p 2 T ), with 1 ≤ m ≤ 2n. This results in a singular behaviour of the cross section as p T → 0. To cure this problem we need to resum these terms to all orders in α S . To properly account for transverse-momentum conservation, the resummation is carried out in impact parameter (b) space [68,69,119]. In this paper we use the formalism of ref. [70]. The partonic transverse-momentum cross section is decomposed as 3) The first term, dσ H,ab , on the right-hand side of eq. (3.3) contains all the logarithmically enhanced terms, and is evaluated by resumming them to all orders. The second term is finite, and can be computed by fixed-order truncation of the perturbative series: it is obtained by computing the standard fixed-order result valid at large p T and subtracting the expansion of the resummed term at the same fixed order. This matching procedure ensures that the resummed and fixed-order components are combined to achieve uniform formal accuracy from the small-to the large-p T region.
The explicit expression of the resummed component is where J 0 (x) is the 0th-order Bessel function and the factor W embodies the all-order resummation of the large logarithmic terms. The all-order structure of W a 1 a 2 is better expressed by considering the N -moments with respect to z = m 2 H /ŝ at fixed m H and is given by . is the Euler number). The function H N in eq. (3.5) does not depend on the impact parameter b and can thus be expanded in powers of α S as JHEP03 (2017)115 where σ 0 (α S , m H ) is the lowest order partonic cross section. The dependence on b is fully contained in the exponent G N (α S ,L, m 2 H /Q 2 ) whose expansion reads where g (1) controls the leading logarithmic (LL) terms, g N the NLL terms, and so forth. The formalism we have briefly recalled defines a systematic expansion of eq. (3.3), whose terms are denoted by NLL+NLO, NNLL+NNLO and so forth. The first label in this notation denotes the logarithmic accuracy, while the second one is referred to the accuracy of the fixed-order calculation. In particular, the NLL+NLO accuracy is obtained by computing the resummed component including the coefficient H (1) together with the functions g (1) and g (2) , and by matching it to the O(α 3 S ) fixed-order result valid at high p T . The NNLL+NNLO accuracy is obtained by including also the coefficient H (2) and the function g  .5), called resummation scale, parametrizes the arbitrariness in the resummation procedure. Its role is analogous to the role played by the renormalization (factorization) scale in the renormalization (factorization) procedure. Although W N does not depend on Q when evaluated at all perturbative orders, an explicit dependence on Q appears when the logarithmic expansion is truncated at a given order. In particular, since the scale Q appears in the definition of the large logarithmic termL, the resummation scale sets the scale up to which the resummation is effective.
As is well known (see section 3 in ref. [70]), the extrapolation of the resummed result at large transverse momenta, where the resummation cannot improve the accuracy of the fixed-order expansion, may lead to unjustified large uncertainties and ensuing lack of predictivity. In the numerical implementation of eq. (3.3) we thus apply a smooth switching procedure as described in ref. [120](see in particular eqs. (13)-(15)). 5 We also point out that, due to the definition of the logarithmic parameter in eq. (3.6), the formalism of ref. [70] fulfils a unitarity constraint such that, upon integration over p T , the customary fixed-order prediction for the inclusive cross section is recovered. More precisely, by performing the resummation at NLL accuracy and including the fixed-order result up to O(α 3 S ) we obtain NLL+NLO accuracy, and the integral of the spectrum is fixed to the NLO total cross section. Despite the switching procedure discussed above the p T spectra we are going to present fulfil such unitary constraint to better than 1%.
Top-and bottom-mass effects can be included in the resummed spectrum along the lines of refs. [72,73]. 6 The inclusion of the top mass does not lead to complications: since m t ∼ m H the computation of the p T spectrum is still a two scale problem. The inclusion of the bottom-mass instead is more difficult. Since m b m H , the computation of the p T spectrum becomes a three scale problem, whose all-order solution is far from being trivial. 7 In ref. [73] a simple solution to this problem was proposed, which implies a choice of different resummation scales for the top and bottom contributions. In particular, since, 5 In the present paper the switching parameters are chosen as p sw T = 125 GeV and ∆pT = 75 GeV. 6 For studies of the resummed pT spectrum in explicit BSM models see for example refs. [121][122][123][124]. 7 For a recent contribution on this subject see ref. [125]. as discussed above, the resummation scale is the scale up to which the resummation is effective, it was suggested to choose for the bottom contribution a lower scale as compared to the top contribution. In ref. [122] this approach has been extended to consider three different resummation scales for the top contribution, the bottom contribution, and the top-bottom interference. We will follow such approach in the next section.
We now discuss the inclusion of BSM effects in the computation. The operators in eqs. (2.4)-(2.6) modify the computation of both the resummed and the fixed-order component in eq. (3.3). However, by limiting ourselves to NLL+NLO, the computation can be greatly simplified. Indeed, the fixed-order result valid at high p T can be obtained by introducing the c t and c b coefficients in the SM amplitude, and supplementing it with an additional contribution, proportional to c g , which corresponds to the QCD amplitude computed in the HTL. As for the resummed component, due to the universality of our formalism [126], its only process dependent contribution is encoded in the coefficients σ 0 and H (1) , which are determined by the Born-level and one-loop amplitudes, respectively. These amplitudes can also be easily obtained from the SM result by introducing the factors c t and c b where appropriate, and adding the point-like HTL amplitude in the SM with a coefficient c g . We emphasize that the first EFT correction to the SM is obtained by interfering the EFT amplitudes with the corresponding SM contributions. With this strategy, we can obtain NLL+NLO predictions for the p T spectrum consistently including the effects of the dimension-six operators. Thanks to the unitarity constraint, the integral of the p T spectrum exactly reproduces the fixed-order NLO result obtained with the inclusion of the same operators.

Results
In this section we present our numerical results for the transverse-momentum spectrum including the effect of dimension-six operators. Our implementation is based on the program HqT [127,128]: a public tool for the computation of the analytic transversemomentum spectrum of the Higgs boson. The contributions from finite top and bottom masses as well as the dimension-six operators are consistently included up to NLL+NLO accuracy. The fixed-order cross section is then cross checked with HIGLU [129] and HNNLO [73,130,131].
We consider pp collisions at √ s = 13 TeV. Our computation is performed in the fiveflavor scheme with the corresponding NLO set of the PDF4LHC2015 [132][133][134][135][136][137] parton distribution functions (PDFs) and the respective value of the strong coupling constant. For the top and bottom quarks running in the loop and for their Yukawa couplings on-shell masses m b = 4.92 GeV and m t = 172.5 GeV are used. As discussed in the previous section, our computation of the NLL+NLO p T spectrum fulfils a unitarity constraint, such that by integrating over p T we recover the fixed-order NLO cross section. The appropriate scale choice for such a resummed computation is of the order of the Higgs boson mass m H . In our study, however, we are also interested in the high-p T region, where a dynamical scale choice has to be preferred. We thus proceed as follows: in the fixed-order computation, which is valid at high p T , our central renormalization and factorization scales are set to In the resummed computation (and its fixed-order expansion) we fix the central scales to µ R = µ F = m H /2. To ensure a proper assignment of the resummation scales for the individual contributions to the cross section we follow the splitting of the SM cross section into a top, a bottom and an interference contribution as suggested in ref. [122] and assign different scales to each of these contributions. In particular we choose: These values are justified by the findings of refs. [73,122,138,139]. As the top contribution is essentially insensitive to the top-quark mass in the small-p T region, where resummation is relevant, we assign Q t also to the contribution with the point-like ggH coupling, when choosing c g = 0. In fact, regarding the splitting of the cross section into the three contributions outlined above we consider the c g amplitude as part of the top amplitude.
In summary, our results for the p T spectrum depend on five scales: the renormalization and factorization scales, and the three resummation scales discussed above. In order to estimate the perturbative uncertainty, we study the corresponding scale variations. As far as renormalization and factorization scales are concerned, the uncertainties are estimated by performing the customary seven-point µ R , µ F variation, i. e. we consider independent variations within the range µ 0 /2 ≤ µ F , µ R ≤ 2µ 0 with 1/2 < µ R /µ F < 2. We then vary each of the three resummation scales (Q t , Q b , Q int ) by a factor of two around their central values (by keeping all the other scales to their central value). We finally combine the ensuing four uncertainty bands by taking the envelope. Figure 2 shows the reference SM prediction together with its perturbative uncertainty. We see that the uncertainty ranges from about ±20% at the peak, to about +50% −30% at p T = 400 GeV.
We start our analysis by considering the individual contribution of exactly one operator. The values of the coefficients c t , c g and c b are varied as much as possible, while requiring the total cross section (integrating over p T ) to not deviate by more than 20% from the SM prediction, which is roughly the current experimental uncertainty on the measured Higgs cross section. 8 Figure 3 shows various predictions of the Higgs transverse-momentum spectrum: SM (black, solid), c t = 1.1 (red, dotted), c t = 0.9 (blue, dashed), c b = 4 (green, dash-dotted), c b = −2 (yellow, short-dashed), c g = 0.008 (magenta, long-dashed) and c g = −0.008 (light-blue, short-dotted). The lower frame illustrates the deviation from the SM prediction by taking the ratio of the curves in the main frame to the black, solid one. The grey-shaded band indicates the uncertainty of the SM result due to scale variations as defined above.
Looking at the low-p T interval (0 GeV≤ p T ≤ 400 GeV) in figure 3 (a) we can directly deduce from the green, dash-dotted and yellow, short-dashed curves that modifications of the bottom Yukawa coupling through c b dominantly affect the low-p T shape of the distribution. In fact, at very low p T we find effects that can even exceed the uncertainty of the SM JHEP03(2017)115 prediction. As expected, c b < 1 (c b > 1) softens (hardens) the spectrum in that region. 9 The point-like Higgs-gluon coupling c g , on the other hand, modifies the p T -shape most notably at large transverse momenta (400 GeV≤ p T ≤ 800 GeV), see figure 3 (b), where a positive (negative) c g value hardens (softens) the spectrum. As expected, modifications of solely the top Yukawa through c t have almost exclusively the effect of a rescaling of the total cross section. By and large all the deviations from the SM prediction through the dimension-six operators are within the scale uncertainty, although the differences in shape give some additional sensitivity to distinguish such effects. It is clear that in order to disentangle effects of this order it is necessary to start from a more accurate SM prediction.
By contrast, the simultaneous variation of more than a single coefficient, as considered in figures 4-6, gives rise to more significant effects. The c g , c t and c b parameters are chosen in the ballpark suggested by the studies of refs. [83][84][85], while still keeping the inclusive cross section within about 20% of its SM value, Indeed, many combinations of c g , c t and c b can be found which mildly affect the total cross section, while significantly changing the shape of the spectrum.
In figure 4 we present the simultaneous variation of c t and c g . The general pattern of these figures follows the pattern of figure 3, but for the variations: c t = 0.1, c g = 0.075 (red, dotted); c t = 0.5, c g = 0.042 (blue, dashed); c t = 1.5, c g = −0.042 (green, dash-dotted) and c t = 2 c g = −0.083 (yellow, short-dashed). In this case, both the small and high-p T     Figure 5. Higgs transverse-momentum spectrum in the SM (black, solid) compared to simultaneous variations of c t and c b for (a) 0 GeV≤ p T ≤ 400 GeV and (b) 400 GeV≤ p T ≤ 800 GeV. The lower frame shows the ratio with respect to the SM prediction. The shaded band in the ratio indicates the uncertainty due to scale variations. See text for more details. behaviour of the spectrum is altered by the different combinations of c t and c g coefficients. It is clear that in particular the large-p T region offers a good discrimination between the different structures of c t and c g in terms of shape. Again, negative (positive) c g values will soften (harden) the spectrum. The effects are well beyond the theoretical uncertainties already at NLL+NLO. We note that the yellow, short-dashed curve corresponding to c t = 2, c g = −0.083 develops a minimum in the ratio to the SM around ∼ 650 GeV. This is due to a compensation between the negative interference between the O 1 and O 2 operators, which is proportional to c g c t and the contribution of O 1 itself, which is proportional to c 2 g and tends to produce a harder spectrum with respect to the SM prediction. Figure 5 shows spectra with modified top and bottom Yukawa couplings: c t = 0.5, c b = −7.46 (red, dotted); c t = 0.8, c b = −3.67 (blue, dashed); c t = 0.9, c b = −1.79 (green, dash-dotted); c t = 1.1 c b = 3.79 (yellow, short-dashed) and c t = 1.2, c b = 4.67 (magenta, long-dashed). In this case, the compensation of the BSM contributions is less straightforward, and it is difficult to compensate c t > 1 without significantly affecting the inclusive cross section. For the magenta, long-dashed curve (c t = 1.2, c b = 4.67) we, thus, allow for a bigger change of the total cross section up to 30%. As pointed out before, the bottom-loop softens the spectrum and, since the variation of the bottom Yukawa coupling is rather large, the squared bottom-term is larger than the top-bottom interference term and the spectrum is softened also when negative c b values are considered. The shape difference to the SM is very significant, but only in the small-p T region, where the softgluon resummation is crucial to obtain a reliable prediction. Indeed, the contribution of the bottom loop decreases with growing p T [140] and above 150 GeV the spectra have all very similar shapes, c t driving their normalization.

JHEP03(2017)115
Finally, we discuss spectra obtained by switching on all three SMEFT operators, as shown in figure 6: c t = 1. Our focus here is on scenarios with increased topquark Yukawa coupling (up to c t = 1.5). These scenarios would be of particular interest in the case in which the excess on the ttH rate over the SM prediction [141,142] should be confirmed. In order to compensate the increase in the cross section driven by c t > 1 a negative c g has been chosen. We observe a general tendency of the BSM spectra to fall below the SM prediction in the intermediate and high transverse-momentum regions, which is due to the negative c g contribution. The total rate is compensated by the enhancement in the low p T region, due to a combination of the negative c g coefficient with both negative and positive c b modifications. Overall, we find sizable distortions of the p T shapes due to the dimension-six operators far beyond the scale uncertainties of the NLL+NLO SM prediction, that exceed the previously considered scenarios with two simultaneous varied coefficients in both size and significance. Despite the similar overall behavior, the predictions for the various scenarios may differ significantly, which enables their discrimination when compared to data.

JHEP03(2017)115
We conclude this section with a comment on the validity of the EFT approach. The computation we have performed is carried out under the assumption that we can consider the effects of higher-dimensional operators as a "small" perturbation with respect to the SM result. This implies in particular that the effect of dimension-eight operators can be neglected. This is not obvious, given that we are studying also the large transversemomentum region. To check the above assumption we have repeated our calculations by dropping the O(1/Λ 4 ) suppressed terms originating from the square of the dimensionsix contributions. We find that in most of the cases the differences with respect to the results shown in figures 3-6 are very small, even at high transverse momenta. Only in the scenarios considered in figure 4 (c t = 0.1, c g = 0.075 and c t = 2, c g = −0.083) the O(1/Λ 4 ) effects are important, and thus, the corresponding quantitative results should be interpreted with care.

Summary
New Physics might be not accessible at the LHC through direct searches, e.g., with the discovery of new resonances. In that case, it is crucial to fully exploit the data to study possible (small) deviations from the SM predictions. SMEFT offers a formalism for the parametrization of high-scale BSM effects, which can be used for this purpose. In the SMEFT framework BSM effects are parametrized through appropriate higher-dimensional operators, and bounds on the corresponding Wilson coefficients can be set by comparing to the experimental data.
In this paper, we have presented a computation of the transverse-momentum spectrum of the Higgs boson in which the SM prediction is supplemented by possible BSM effects. Such effects are modeled by augmenting the SM Lagrangian with appropriate dimensionsix operators. Our calculation consistently includes all the terms up to O(α 3 S ) accuracy and is supplemented by soft-gluon resummation at NLL accuracy, which is required to obtain reliable predictions at small transverse momenta. At the same level of accuracy we implement three dimension-six operators, related to the modifications of top and bottom Yukawa couplings and to the inclusion of a point-like ggH coupling. Additionally, we studied the impact of the chromomagnetic operator on the Higgs cross section at LO, which had been previously addressed in the literature by different groups with contradicting results.
We have constructed a tool for reliable predictions of the Higgs p T distribution including dimension-six operators and performed a comprehensive study of the possible effects due to the different dimension-six operators, by studying the impact of variations of c t , c b and c g on the transverse-momentum spectrum of the Higgs boson. We varied the above coefficients in the range suggested by recent global analyses and required the total cross section to meet the SM prediction at NLO within the current O(20%) experimental uncertainty. Our results can be summarized as follows: coupling of the Higgs boson to gluons (O 1 ), on the other hand, changes the shape of the distribution most notably in the high-p T tail. As expected, changes in the topquark Yukawa coupling (O 2 ) primarily affect the normalization and approximately correspond to a simple rescaling of the spectrum.
• The shape of the transverse momentum distribution depends on the mass of the particle that mediates the Higgs-gluon coupling. The lower the mass of that particle, the softer is the resulting spectrum. Therefore, the p T shape associated with the bottom loop is softer, in particular at small transverse momenta, than the SM one and, when increasing the absolute value of the bottom-quark Yukawa coupling, positive (negative) values soften (harden) the spectrum, if the top-bottom interference is dominant (small variations of c b ). In contrast the spectrum becomes always softer for |c b | 1, if the squared bottom-loop is dominant (large variations of c b ). Furthermore, a point-like coupling between gluons and the Higgs boson leads to the hardest spectrum and a positive (negative) c g value hardens (softens) the shape as compared to a Higgs boson mediated by a top-quark loop.
• While individual variations of the various Wilson coefficients produce effects that hardly exceed the NLL+NLO perturbative uncertainties, the simultaneous variation of two or more operators can significantly distort the spectrum, still keeping the total rate consistent with the NLO prediction within the current experimental uncertainties.
• Choosing combinations of c t , c b and c g that compensate each other at the level of the total cross section allows us to deform the shape of the Higgs p T spectrum far beyond the uncertainties of our NLL+NLO prediction in the SM and thus allow for a disentanglement of the different Wilson coefficients in future analyses at the LHC.
We conclude our discussion by adding a few comments on the limitations of the calculation presented here. When only small deviations from the SM are considered the theoretical uncertainty affecting the SM prediction becomes relevant. We have seen that at NLL+NLO the uncertainties from missing higher-order contributions, estimated through scale variations, are about ±20% at the peak and increase by roughly a factor of two at high transverse momenta. These uncertainties imply an unavoidable limitation for the extraction of constraints on the Wilson coefficients and it is only through their reduction that the sensitivity to BSM effects can be increased. The natural question is whether the calculation we have carried out can be extended to the next order, i.e., to NNLL+NNLO. To consistently carry out such extension we would need the heavy-quark mass effects at NNLO, which are currently unavailable. A possible way out is to include the effects beyond NLL+NLO in the HTL, as is currently done in state of the art SM calculations [73]. This approach implies that the relevant Higgs production amplitudes would contain a c g term already in the SM. Nonetheless, it is reasonable to assume that the QCD effects beyond our NLL+NLO accuracy factorise with respect to the BSM corrections. In this approximation, the relative BSM/SM effects we have obtained in this paper (i.e., the ratios plotted -15 -

JHEP03(2017)115
in the lower panels of figures 3-6) can be directly used to include BSM effects on top of NNLL+NNLO accurate SM predictions.
Another aspect which deserves some comments is the set of dimension-six operators we have considered. In the present calculation we have limited ourselves to consider the contributions of the operators related to modified top and bottom Yukawa couplings and of the additional direct Hgg interaction. As discussed in section 2, although formally suppressed by two powers of the top Yukawa coupling, the chromomagnetic operator could still significantly contribute, within the current bounds, to the gluon fusion cross section. The extension of our calculation to include these effects is left for future work.