# Matching uncertainties in the prediction of the Higgs boson transverse
momentum in the SM and beyond

Deutsches Elektronen-Synchrotron (DESY), Notkestraße 85, D-22607 Hamburg,
Germany
E-mail

###### Abstract:

We present the results of our recent study [1] of the theoretical
uncertainties that affect the predictions for the Higgs-boson transverse-
momentum in gluon fusion when fixed- and all-order results are matched. Our
investigation consists of a twofold analysis: first we present a detailed
comparison of two recently introduced prescriptions for the determination of
the matching scale [2, 3], then we apply the results of these methods to three
widely used matching frameworks, namely the aMC@NLO and POWHEG Monte Carlo
approaches and analytic resummation. The results of our study are applied to
the production of the SM Higgs boson and of the neutral Higgs bosons of the
Two-Higgs-Doublet Model in a variety of scenarios.

††dedicated: DESY 16-173

##  1 Introduction

After the discovery of the Higgs boson at the LHC [4, 5], a quest has started
to characterize the newly found particle. Its properties have already been a
topic of studies of many articles whose aim was to understand the
compatibility with the Standard Model (SM) predictions. So far these studies
have relied on the precise predictions available for production and decay
rates. Only in the last year the experimental measurements for differential
observables have been published [6, 7]. As already stressed by various authors
[8, 9, 10, 11, 12, 13, 14, 15, 16] the Higgs-transverse momentum ($p_{\bot}$)
distribution opens the possibility of probing the loop dynamic of the gluon
fusion process. This observable is therefore sensible to modifications of the
Yukawa couplings and/or to the presence of new states beyond the SM ones. The
opportunity of exploiting the $p_{\bot}$ measurement is becoming more and more
interesting as the LHC accumulates a new wealth of data during the second run
of operations. At the same time, when considering models of new physics with
enlarged Higgs sectors, an accurate description of the transverse momentum,
which can show sensible deviations from the SM prediction, is required.

To properly describe the transverse momentum distribution, whose fixed-order
prediction is logarithmically divergent in the limit $p_{\bot}\to 0$, one
needs to resum the terms enhanced by powers of $\log(p_{\bot}/m_{h})$ to all
orders in $\alpha_{s}$, where $m_{h}$ is the mass of the scalar resonance.
This resummation is usually performed either analytically or algorithmically.
Due to the theoretical formalism on which it is based, the resummation
procedure is strictly valid only in the limit of collinear emissions and
therefore, for gluon fusion, in the limit of zero transverse momentum of the
Higgs boson. Therefore, to properly describe the whole $p_{\bot}$-spectrum, a
matching between the fixed order and the resummed results is needed.
Particular care has to be taken to avoid any kind of double counting. This has
been achieved in various frameworks, both analytic [17, 18, 19, 20] and
numeric [21, 22]. Common to all the approaches is the introduction of a new
unphysical scale, which we will subsequently denote as the _matching scale_
($\mu$), whose role is to define the transition region between the fixed- and
the all-order results. The dependence of the matched result on this scale is
of higher logarithmic order, however a careless choice can ruin the
perturbative convergence of the result. This is especially true for those
processes that are characterized by more than one scale, as the one that we
are considering in our study, as we will see in the next sections.

##  2 The Higgs transverse momentum as a multiscale problem

Higgs production in gluon fusion has been originally studied in the so-called
Heavy Quark Effective Field Theory (HQEFT) obtained in the limit where the
top-quark mass is taken to be very large compared to Higgs boson mass. The use
of the HQEFT has the advantage of reducing what is a one-loop LO process to a
tree level one, with a sensible decrease of the complexity of the computation.
Under this approximation, the total cross section has been computed up to
$\mathrm{N}^{3}\mathrm{LO}$, while differential computations are available up
to NNLO (see ref. [23, 24, 25] for a complete reference of all the results
available). However, being this an effective description, it is valid only if
we are not probing mass scales that are equal or larger than the top quark
mass. Therefore, the HQEFT is a description limited to Higgs boson masses
smaller than the top quark mass, for what concerns the total cross section,
and to $p_{\bot}$ less than top mass for the transverse momentum distribution.
Moreover, it neglects completely the contribution coming from diagrams where
the coupling of the Higgs to the gluons is mediated by a loop of light quarks.
The latter are important for precision predictions in the SM and fundamental
for Beyond-Standard-Model (BSM) Higgs boson production, where it can happen
for specific regions of the BSM parameter space that the dominant contribution
to the cross section is from bottom quark diagrams, oppositely to the SM.
Therefore, to properly describe the Higgs boson transverse momentum, we need
to perform the computation in the full theory, being it either the SM or a BSM
model as the THDM or the MSSM. In the case of the SM and the THDM, complete
computations are available up to NLO. In the case of the MSSM, even the full
NLO result is not known analytically.

Restricting ourselves, for simplicity to a THDM-like scenario, the description
of the Higgs transverse momentum in gluon fusion is characterized by four mass
scales: the Higgs boson mass; the top quark mass; the bottom quark mass; the
$p_{\bot}$ of the radiated parton. All these physical scales and their non-
trivial interplay have to be taken into account in the choice of the matching
scale.

To follow this requirement, it was first proposed by Grazzini et al. [26] to
split the complete squared matrix element into components that are
characterized just by a single mass scale or by a specific combination
thereof. Originally this split was into two parts, the top contribution and
the bottom plus the interference terms. Here we follow more recent
developments where the amplitude is divided into three terms [2, 3]. To
achieve this, we rewrite the full amplitude as

$\displaystyle|\mathcal{M}(t+b)|^{2}=\overbrace{|\mathcal{M}(t)|^{2}}^{\mathrm{only~{}top}}+\overbrace{|\mathcal{M}(b)|^{2}}^{\mathrm{only~{}bottom}}+\underbrace{\left(|\mathcal{M}(t+b)|^{2}-|\mathcal{M}(t)|^{2}-|\mathcal{M}(b)|^{2}\right)}_{\mathrm{interference}}\,,$
(1)

then compute separately the matched prediction for each contribution, with a
different matching scale each, and finally sum all the three components
together. We note that, while the decomposition given in eq. (1) is a trivial
identity at the level of the total cross section, because the latter is
independent on the matching scale, it yields a modified shape for the
transverse momentum distribution with respect to the one obtained using a
single scale for the full amplitude. Our master formula for the best
prediction of the Higgs boson transverse momentum distribution is therefore
given by

$\displaystyle\frac{d\sigma}{dp_{\bot}}=\frac{d\sigma_{t}}{dp_{\bot}}\bigg{|}_{\mu_{t}}+\frac{d\sigma_{b}}{dp_{\bot}}\bigg{|}_{\mu_{b}}+\frac{d\sigma_{int}}{dp_{\bot}}\bigg{|}_{\mu_{\mathrm{int}}}\,,$
(2)

where with $\mu_{t}$, $\mu_{b}$ and $\mu_{\mathrm{int}}$ we have denoted
respectively the matching scale for the top and bottom contributions and the
interference term.

###  2.1 Matching scale determination

Figure 1: On the left, ratio ($C$) between the full squared matrix amplitude
and the collinear limit of the gluon initiated subprocess for the top, bottom
and interference terms as well as for the HQEFT computation, as a function of
the Higgs transverse momentum; the dashed lines indicate a deviation of $\pm
10\%$ from one. On the right, ratio between the matched and the fixed order
predictions at the hadronic level, for the top contribution and for different
values of matching scale. The latter is here denoted by
$\mathrm{Q}_{\mathrm{res}}$ because the results shown are obtained in the
framework of AR.

In our work we have compared two prescriptions that recently appeared in the
literature111In the context of SCET, though only in the HQEFT, a detailed
study on the problem of scale determination has been published in ref. [30,
31, 32].: the partonic analysis published in ref. [3] (BV) and the hadronic
analysis in ref. [2] (HMW).

The BV prescription is based upon the observation that the resummation
formalism relies on the factorization of the squared matrix element in the
limit of soft and/or collinear emissions (see also ref. [27]). In detail, the
accuracy of the collinear approximation of the gluon fusion process, as a
function of $p_{\bot}$, is evaluated at the partonic level, separately for
each term in eq. (1). If the approximation is violated by a well-defined
threshold, here chosen to be $10\%$, for a given value of $p_{\bot}$, then the
latter is taken to be the matching scale to be used in the matched computation
(see the left plot of fig. 1). The procedure is applied separately to the
gluon-gluon and the quark-gluon subprocesses. The results are then averaged
with a differential-weight that keeps into account that the two channels
contribute with varying proportions at different $p_{\bot}$ due to the
distinct behavior of the quark and gluon PDFs.

The HMW method follows from two principles: for $p_{\bot}\gtrsim m_{h}$ the
spectrum is correctly described by fixed-order perturbation theory and
therefore the latter should be the definitive prediction in this range; one
would like to have an all-order result for a $p_{\bot}$-range as large as
possible. Practically, this translates into the definition of a scale
$Q^{max}_{res}$ as the maximum scale for which the resummed distribution is
within the interval $[0,2]\cdot[d\sigma^{fNLO}/dp_{\bot}]$ for $p_{\bot}\geq
m_{h}$. The matching scale $\mu$ is then taken to be equal to half
$Q^{max}_{res}$ (see the right plot of fig. 1).

We stress that, in both cases, the scales are independent on the Yukawa
coupling of the quark to the Higgs.

###  2.2 Matching scale comparison

Figure 2: Comparison of the matching scales, for the top (red) and bottom
(green) contributions and for the interference term (blue), obtained with
either the BV ($w_{i}$, dashed lines) or the HMW prescription ($Q_{i}$, solid
lines) as a function of the Higgs mass.

In fig. 2 we show the two scale sets as a function of the Higgs boson mass,
for the scalar case. Due to the different assumptions on which the two
procedures are based, it is not surprising that the numerical values are
different. Indeed the BV prescription is sensitive to the behavior of the
transverse spectrum in the low-$p_{\bot}$ region, while the HMW method is
designed around the high-$p_{\bot}$ tail. What we observe is a moderate
agreement for the top contribution (with only the BV scale showing a
sensitivity to the $2\cdot m_{\mathrm{top}}$ threshold) and a very good
agreement for the bottom one. Instead, for the interference term the two
prescriptions can differ by a large amount and specifically this happens when
the LO term is much smaller than the NLO one. Indeed, in this case we have
that the resummed contribution, being proportional222Apart from corrections
due to the virtual contributions that are small compared to the total cross
section. to the LO term, is small and the Higgs transverse momentum will be
given almost completely by the hard emission from the NLO term. Then the
collinear approximation will fail for any value of $p_{\bot}>0$ and therefore
the BV scale will vanish; on the contrary, because the matched curve is almost
identical to the fixed order one, for every value of the resummation scale,
the HMW scale will tend to be very large. As a general feature of both scale
sets, we observe that, for heavy Higgs masses, the scales for all the three
contributions are much smaller than the commonly used value of $m_{h}/2$.

##  3 Simulation setup

In our study of the theoretical uncertainties of the Higgs transverse momentum
distribution we used three different codes: MoRe-SusHi, which implements the
Analytic Resummation (AR) procedure; gg_H_2HDM from the POWHEG-BOX framework;
aMCSusHi, based upon the aMC@NLO framework and the SusHi [28] amplitudes. The
uncertainty band due the matching is determined using the following
prescription: given the reference values
$(\mu_{t},\mu_{b},\mu_{\mathrm{int}})$ for the three matching scales discussed
in the previous sections, we consider all the possible combinations generated
by taking half and twice these values or the reference values themselves; for
each combination we compute our prediction for the Higgs $p_{\bot}$; finally
we take the envelope for each $p_{\bot}$ bin, i.e. we take the maximum and the
minimum value among all the predictions. Only for AR, we follow ref. [2] and
we apply an additional damping factor to the error band at large $p_{\bot}$.

In our study, besides the SM, we considered various THDM-II scenarios. Each of
the latter was chosen because it is characterized by a specific feature, e.g.
a very large Yukawa coupling to one of the two quarks, whose impact on the
$p_{\bot}$ distribution we want to understand. We also considered both scalar
and pseudoscalar productions. In this proceeding, due to the restricted space
available, we report the results only for the SM Higgs boson using BV scales
and for heavy scalar production in the large-$b$ scenario with the HMW scale
set. See table 1 of ref. [1] for a list of all scenarios.

All the numerical results are computed for the LHC, with a a center-of-mass
energy of $\sqrt{S}=13$ TeV, using the MSTW2008nlo68cl PDF set and the
associated value of $\alpha_{s}(M_{Z})=0.120179$. The renormalization and
factorization scales are set to the Higgs boson mass; the quark masses are
fixed at $m_{\mathrm{top}}=172.5$ GeV and $m_{b}=4.75$ GeV respectively; the
Yukawa couplings of the Higgs to quarks are renormalized in the On-Shell (OS)
scheme. We used Pythia8 as our Parton Shower [29].

##  4 SM results

Figure 3: Shapes of the transverse momentum distribution for a SM Higgs boson
of mass $m_{h}=125$ GeV. We show the curves from AR (black, solid), MC@NLO
(red, dotted) and POWHEG (blue, dashed overlaid by points) for the BV scale
set. For reference we also show the fixed order result (green, dash-dotted
with open boxes). On the left, we plot the spectrum up to $400$ GeV, on the
right we show a zoom around the first $100$ GeV. In the main frame we show the
absolute distributions, in the middle inset the ratio of the central curves to
the analytic resummation result and in the bottom inset the uncertainty bands,
again normalized to the AR resummation value.

In fig. 3 we show the shape of the transverse momentum distribution (i.e. the
integral of each curve is normalized to one) obtained using the BV scale set
for a SM Higgs boson of $m_{h}=125$ GeV. Qualitatively we observe that all the
codes agree within their uncertainty bands, at least for $p_{\bot}<m_{h}$.
More in detail, we observe that the two MC event generators, POWHEG and MC@NLO
are in excellent agreement in the region $10<p_{\bot}/\mathrm{GeV}<130$, while
they differ by about $20\%$ from the central AR prediction. The position of
the peak is also slightly different between the MCs and AR. Turning now our
attention to the high-$p_{\bot}$ tail, we see that: the AR prediction
approaches the fNLO at the level of $5\%$ above $p_{\bot}\simeq 130$ GeV; the
transition to the NLO prediction in MC@NLO is around $p_{\bot}\simeq 180$ GeV;
POWHEG, on the other hand, always remains $20\%$ above the fixed order result.
The latter is a general feature of the POWHEG matching that we observe and
that will be analyzed more in detail in the THDM analysis that follows.

Concerning the uncertainty bands, for $p_{\bot}>130$ GeV AR has no uncertainty
band due to the artificial suppression given by the damping factor, as
mentioned in section 3; on the other hand, the uncertainty for the MC@NLO
prediction is of order $\pm 10\%$, while the width of the POWHEG band
decreases uniformly from $\pm 20\%$ to $\pm 10\%$. In the intermediate region,
all three codes show a bulgy structure with a maximum of $\pm 20\%$ and $\pm
35\%$ for AR and the MCs respectively and a minimum of just a few percent
above the peak position. Finally, in the small-$p_{\bot}$ region, the AR
uncertainty band grows to $100\%$, the POWHEG one to $40\%$, while MC@NLO
shows only a moderate $\pm 15\%$ uncertainty.

##  5 THDM results

Figure 4: The plot on the left is as in fig. 3 but for the production of a
heavy CP-even Higgs of $m_{H}=300$ GeV and for the use of the HMW scales. On
the right we show the results with the MC@NLO and POWHEG curves obtained with
the changes to the shower scale prescription as described in the text.

The left plot of fig. 4 shows the same curves as the SM case discussed in the
previous section. Differently from the SM, however, the discrepancies between
the three frameworks are more marked. Concerning the central curves, we
observe that POWHEG produces a spectrum that is significantly harder, being
over $50\%$ above the fNLO result for $p_{\bot}>200$ GeV. On the other hand,
in the intermediate region between $10\lesssim p_{\bot}/\mathrm{GeV}\lesssim
130$, POWHEG and AR agree within $10\%$ while the MC@NLO curve is
substantially larger. At smaller transverse momentum, for $p_{\bot}<30$ GeV,
the two MCs have a much better agreement.

The behavior of the uncertainty bands is also quite different in the various
frameworks: the MC@NLO band blows up to $\mathcal{O}(100\%)$ around
$p_{\bot}\simeq 125$ GeV; the POWHEG band remains very small all over the
whole $p_{\bot}$ range.

To understand the origin of these differences, we investigated the dependence
of the MC predictions on the prescription for the shower scale $Q_{sh}$. The
latter is the scale that it is passed to the PS to be used as an upper bound
for the $p_{\bot}$ of the emitted radiation. Concerning POWHEG, we notice that
by restricting, for the class of events333These are the _remnant_ events. By
default setting of the POWHEG-BOX their shower scale is set to the transverse
momentum of the emitted parton. See ref. [33] for a detailed description of
the POWHEG-BOX. that describe the high-$p_{\bot}$ tail, the shower scale to be
at most the matching scale (either BV or HMW), we recover the fixed order
result in the same way as AR and MC@NLO do (as it can be seen from the purple
curve on the right plot of fig. 4). Moreover, the shape of the uncertainty
band is also changed, showing now a bulge between $50$ GeV and $100$ GeV.
Relatively to MC@NLO, we first recall that by default the shower scale is
extracted from a probability distribution dependent on the LO kinematic and
centered around the matching scale [34]; if we replace the default
distribution with the $\delta$-function $\delta(Q_{sh}-\mu_{i})$, we observe a
significant change for both the central prediction and uncertainty bands, as
it can be seen by the yellow curve in the right plot of fig. 4. These
observations lead us to conclude that for this specific scenario there is a
high sensitivity not only to the numerical values of the scales but also to
the specific details of the matching procedure.

##  6 Conclusions

In this talk we have presented the results of our recent study [1] of the
theoretical uncertainties intrinsic to the matching procedure between fixed-
and all-order results in the computation of the transverse momentum
distribution of the Higgs boson in gluon fusion. Specifically for this
process, which involves different mass scales, even the choice of the central
values for the matching scales has become a matter of debate. In this context,
we performed a thorough analysis of the predictions obtained using three
different matching frameworks (analytic resummation, POWHEG and MC@NLO) and
two different prescriptions for the determination of the matching scales (BV
or HMW). Our comparison was twofold: first we addressed the issue of the
determination of the central value for the matching scale for the top and
bottom contributions and for the interference term, by providing a qualitative
and quantitative comparison of the BV and HMW approaches; then we compared the
results for the shape of the $p_{\bot}$ distribution obtained with different
scale-choices and frameworks.

We have found that the prediction of the Higgs transverse momentum is affected
by uncertainties444In this study we have limited ourselves to the matching
uncertainty, however the latter should always be combined with the fixed-order
perturbative uncertainty, usually estimated through the variation of the
renormalization and factorization scale. up to several tens of percent,
depending on the scenario, the $p_{\bot}$ value and the framework under
consideration. In the low-$p_{\bot}$ region we find reasonable agreement
between the different codes, although AR usually shows a much softer spectrum.
However, in all the three frameworks and especially for AR, the error bands
grow in this region, therefore providing compatibility between the different
results. In the intermediate region, we find non-trivial differences between
the three frameworks, which are more pronounced in the bottom dominated
scenarios. In the latter case, we also find a large dependence on the specific
details of the matching formulation inside each framework. In the
large-$p_{\bot}$ tail, where technically all the codes have LO accuracy 555In
the SM case, new developments are available which provide NLO-QCD accuracy in
the description of the Higgs $p_{\bot}$, see refs. [35, 36, 37, 38, 39]., we
find that POWHEG systematically predicts a harder spectrum than MC@NLO and AR.
The latter are instead softer and compatible with the fNLO result. We
identified one source of this difference in the prescription used to define
the allowed phase space for radiation emission by the PS. Restricting the
phase space, as it happens already in the MC@NLO framework, allows also POWHEG
to approach the fNLO at high-$p_{\bot}$, as we have shown with a dedicated
analysis.

## References

  * [1] E. Bagnaschi, R. V. Harlander, H. Mantler, A. Vicini and M. Wiesemann, JHEP 1601 (2016) 090 doi:10.1007/JHEP01(2016)090 [arXiv:1510.08850 [hep-ph]].
  * [2] R. V. Harlander, H. Mantler and M. Wiesemann, JHEP 1411 (2014) 116 doi:10.1007/JHEP11(2014)116 [arXiv:1409.0531 [hep-ph]].
  * [3] E. Bagnaschi and A. Vicini, JHEP 1601 (2016) 056 doi:10.1007/JHEP01(2016)056 [arXiv:1505.00735 [hep-ph]].
  * [4] G. Aad et al. [ATLAS Collaboration], Phys. Lett. B 716 (2012) 1 doi:10.1016/j.physletb.2012.08.020 [arXiv:1207.7214 [hep-ex]].
  * [5] S. Chatrchyan et al. [CMS Collaboration], Phys. Lett. B 716 (2012) 30 doi:10.1016/j.physletb.2012.08.021 [arXiv:1207.7235 [hep-ex]].
  * [6] G. Aad et al. [ATLAS Collaboration], Phys. Rev. Lett. 115 (2015) no.9, 091801 doi:10.1103/PhysRevLett.115.091801 [arXiv:1504.05833 [hep-ex]].
  * [7] V. Khachatryan et al. [CMS Collaboration], Eur. Phys. J. C 76, no. 1, 13 (2016) doi:10.1140/epjc/s10052-015-3853-3 [arXiv:1508.07819 [hep-ex]].
  * [8] U. Langenegger, M. Spira, A. Starodumov and P. Trueb, JHEP 0606 (2006) 035 doi:10.1088/1126-6708/2006/06/035 [hep-ph/0604156].
  * [9] O. Brein and W. Hollik, Phys. Rev. D 76 (2007) 035002 doi:10.1103/PhysRevD.76.035002 [arXiv:0705.2744 [hep-ph]].
  * [10] E. Bagnaschi, G. Degrassi, P. Slavich and A. Vicini, JHEP 1202 (2012) 088 doi:10.1007/JHEP02(2012)088 [arXiv:1111.2854 [hep-ph]].
  * [11] R. V. Harlander and T. Neumann, Phys. Rev. D 88 (2013) 074015 doi:10.1103/PhysRevD.88.074015 [arXiv:1308.2225 [hep-ph]].
  * [12] C. Grojean, E. Salvioni, M. Schlaffer and A. Weiler, JHEP 1405 (2014) 022 doi:10.1007/JHEP05(2014)022 [arXiv:1312.3317 [hep-ph]].
  * [13] A. Azatov and A. Paul, JHEP 1401 (2014) 014 doi:10.1007/JHEP01(2014)014 [arXiv:1309.5273 [hep-ph]].
  * [14] A. Banfi, A. Martin and V. Sanz, JHEP 1408 (2014) 053 doi:10.1007/JHEP08(2014)053 [arXiv:1308.4771 [hep-ph]].
  * [15] S. Dawson, I. M. Lewis and M. Zeng, Phys. Rev. D 90 (2014) no.9, 093007 doi:10.1103/PhysRevD.90.093007 [arXiv:1409.6299 [hep-ph]].
  * [16] U. Langenegger, M. Spira and I. Strebel, arXiv:1507.01373 [hep-ph].
  * [17] J. C. Collins, D. E. Soper and G. F. Sterman, Nucl. Phys. B 250 (1985) 199. doi:10.1016/0550-3213(85)90479-1
  * [18] G. Bozzi, S. Catani, D. de Florian and M. Grazzini, Nucl. Phys. B 737 (2006) 73 doi:10.1016/j.nuclphysb.2005.12.022 [hep-ph/0508068].
  * [19] S. Mantry and F. Petriello, Phys. Rev. D 81 (2010) 093007 doi:10.1103/PhysRevD.81.093007 [arXiv:0911.4135 [hep-ph]].
  * [20] T. Becher, M. Neubert and D. Wilhelm, JHEP 1305 (2013) 110 doi:10.1007/JHEP05(2013)110 [arXiv:1212.2621 [hep-ph]].
  * [21] S. Frixione and B. R. Webber, JHEP 0206 (2002) 029 doi:10.1088/1126-6708/2002/06/029 [hep-ph/0204244].
  * [22] P. Nason, JHEP 0411, 040 (2004) doi:10.1088/1126-6708/2004/11/040 [hep-ph/0409146].
  * [23] S. Dittmaier et al. [LHC Higgs Cross Section Working Group Collaboration], doi:10.5170/CERN-2011-002 arXiv:1101.0593 [hep-ph].
  * [24] S. Dittmaier et al., doi:10.5170/CERN-2012-002 arXiv:1201.3084 [hep-ph].
  * [25] S. Heinemeyer et al. [LHC Higgs Cross Section Working Group Collaboration], doi:10.5170/CERN-2013-004 arXiv:1307.1347 [hep-ph].
  * [26] M. Grazzini and H. Sargsyan, JHEP 1309 (2013) 129 doi:10.1007/JHEP09(2013)129 [arXiv:1306.4581 [hep-ph]].
  * [27] A. Banfi, P. F. Monni and G. Zanderighi, JHEP 1401 (2014) 097 doi:10.1007/JHEP01(2014)097 [arXiv:1308.4634 [hep-ph]].
  * [28] R. V. Harlander, S. Liebler and H. Mantler, Comput. Phys. Commun. 184 (2013) 1605 doi:10.1016/j.cpc.2013.02.006 [arXiv:1212.3249 [hep-ph]].
  * [29] T. Sjostrand, S. Mrenna and P. Z. Skands, Comput. Phys. Commun. 178 (2008) 852 doi:10.1016/j.cpc.2008.01.036 [arXiv:0710.3820 [hep-ph]].
  * [30] Z. Ligeti, I. W. Stewart and F. J. Tackmann, Phys. Rev. D 78 (2008) 114014 doi:10.1103/PhysRevD.78.114014 [arXiv:0807.1926 [hep-ph]].
  * [31] R. Abbate, M. Fickinger, A. H. Hoang, V. Mateu and I. W. Stewart, Phys. Rev. D 83 (2011) 074021 doi:10.1103/PhysRevD.83.074021 [arXiv:1006.3080 [hep-ph]].
  * [32] C. F. Berger, C. Marcantonini, I. W. Stewart, F. J. Tackmann and W. J. Waalewijn, JHEP 1104 (2011) 092 doi:10.1007/JHEP04(2011)092 [arXiv:1012.4480 [hep-ph]].
  * [33] S. Alioli, P. Nason, C. Oleari and E. Re, JHEP 1006 (2010) 043 doi:10.1007/JHEP06(2010)043 [arXiv:1002.2581 [hep-ph]].
  * [34] J. Alwall et al., JHEP 1407 (2014) 079 doi:10.1007/JHEP07(2014)079 [arXiv:1405.0301 [hep-ph]].
  * [35] K. Hamilton, P. Nason, E. Re and G. Zanderighi, JHEP 1310 (2013) 222 doi:10.1007/JHEP10(2013)222 [arXiv:1309.0017 [hep-ph]].
  * [36] S. Hoeche, Y. Li and S. Prestel, Phys. Rev. D 90 (2014) no.5, 054011 doi:10.1103/PhysRevD.90.054011 [arXiv:1407.3773 [hep-ph]].
  * [37] K. Hamilton, P. Nason and G. Zanderighi, JHEP 1505 (2015) 140 doi:10.1007/JHEP05(2015)140 [arXiv:1501.04637 [hep-ph]].
  * [38] S. Alioli, C. W. Bauer, C. Berggren, F. J. Tackmann, J. R. Walsh and S. Zuberi, JHEP 1406 (2014) 089 doi:10.1007/JHEP06(2014)089 [arXiv:1311.0286 [hep-ph]].
  * [39] R. Frederix, S. Frixione, E. Vryonidou and M. Wiesemann, JHEP 1608, 006 (2016) doi:10.1007/JHEP08(2016)006 [arXiv:1604.03017 [hep-ph]].