Refining the GENEVA method for Higgs boson production via gluon fusion

We describe a number of improvements to the GENEVA method for matching NNLO calculations to parton shower programs. In particular, we detail changes to the resummed calculation used in the matching procedure, including disentangling the cross section dependence on factorisation and beam scales, and an improved treatment of timelike logarithms. We also discuss modifications in the implementation of the splitting functions which serve to make the resummed calculation differential in the higher multiplicity phase space. These changes improve the stability of the numerical cancellation of the nonsingular term at small values of the resolution parameter. As a case study, we consider the gluon-initiated Higgs boson production process $gg\to H$. We validate the NNLO accuracy of our predictions against independent calculations, and compare our showered and hadronised results with recent data taken at the ATLAS and CMS experiments in the diphoton decay channel, finding good agreement.


Introduction
In recent years, the quest for precision at the Large Hadron Collider (LHC) has seen many impressive milestones in the development of theoretical tools used to describe hadronic collisions. Many processes are currently known at next-to-next-to-leading order (NNLO) in perturbative QCD, and several 2 → 1 processes even at one order higher (N 3 LO). One particular direction in which much fruitful progress has been made is in the matching of higher order perturbative calculations to parton shower (PS) programs, resulting in Monte Carlo event generators which combine the advantages of fixed-order calculations with the flexibility of parton shower tools. This paradigm generally goes under the name of NNLO+PS.
Several different methods which reach NNLO+PS accuracy have been proposed [1][2][3][4][5][6][7], with most applications using a resummed calculation -either directly or via the Sudakov factor in a shower Monte Carlo -in a suitable resolution variable alongside the fixed order to achieve the matching. Of these methods, the Geneva approach [4,8] has the advantage of being particularly flexible with regard to the framework used for the resummed calculation and the choice of the resolution variable, while also exploiting the possibility of reaching higher logarithmic accuracies in both direct QCD and soft-collinear effective theory (SCET) formalisms. This has resulted in the application of the method to a number of colour singlet processes [4,[9][10][11][12][13][14][15], as well as first steps towards implementations involving coloured final states [16].
In this work, we describe a number of improvements to the Geneva event generator, which both extend the capabilities of the program and improve its numerical performance. Specifically, we detail a new treatment of the splitting functions, which were first introduced in the original Geneva implementation [8] and serve to make the resummed calculation used in the matching procedure differential in the higher multiplicity phase space. The new approach significantly increases the performance of the code in extreme soft and collinear regions, where the cancellation of large logarithmic terms is extremely delicate. We also implement a more rigorous treatment of the theoretical uncertainties by disentangling the factorisation and renormalisation scale dependences in the cross section and allowing their independent variations. This puts our uncertainty estimation on a robust and more conservative theoretical footing, and will also prove important for the implementation of processes featuring perturbatively generated heavy flavours in the initial state. Finally, we also discuss the treatment of timelike logarithms in our calculation, the inclusion of which has been shown to improve the perturbative convergence of colour singlet production processes such as gg → H [17][18][19].
In order to study the various improvements to the Geneva program, we have implemented the gluon-initiated Higgs boson production process (gg → H). We use a resummed calculation in the zero-jettiness resolution variable T 0 obtained via SCET up to N 3 LL accuracy. The process is interesting from both an experimental and a theoretical perspective for a number of reasons.
From the theory side, many calculations work in the limit in which the top-quark mass is considered to be large compared to other scales present in the process, the so-called heavy-top limit (HTL). This significantly simplifies the computational complexity since the top-quark loop coupling the Higgs boson to gluons is integrated out, resulting in an effective ggH vertex and further effective vertices with more gluons and Higgs bosons. Consequently, calculations including QCD corrections up to N 3 LO are now available in this limit [36][37][38][39][40][41][42][43], including matching to resummed calculations up to N 3 LL ′ accuracy in transverse momentum [43][44][45][46][47][48][49][50][51][52] and in jet veto observables [53][54][55][56][57]. There has also been a considerable amount of work on improving calculations beyond the HTL by including quark mass effects [58][59][60][61][62]; this has culminated in a calculation of the exact top-quark mass dependence at NNLO in QCD [63,64]. Additionally, the fact that the perturbative series is known to be poorly convergent has motivated the study of alternative scale choices which include π 2 terms arising from kinematic logarithms at all orders. Finally, the simplicity of this process in terms of its kinematics and matrix elements makes it a particularly suitable testing ground for the improvements which we will detail in this work.
The rest of the paper is organised as follows. In sec. 2, we provide a brief recap of the Geneva method and its application to gluon-induced Higgs production, before discussing the new features which have been implemented in the program in sec. 3. In sec. 4 we validate the NNLO accuracy of our calculation for the gg → H process and discuss the matching to the parton shower provided by Pythia8 [65]. Finally, we compare our results with the pp → H → γγ data collected at the ATLAS and CMS experiments [66,67] in sec. 5. We give our conclusions in sec. 6. In App. A we show an application of the novel splitting function implementation to the Drell-Yan process.

Theoretical framework
In the following we lay out the theoretical framework we work in. We start by giving a summary of the Geneva event generator formalism, which includes the matching procedure of the fixed-order calculation to the resummed prediction. We then focus on the definition of the process under study, i.e. Higgs boson production via gluon fusion, and on its zerojettiness resummation.

The GENEVA method
The complete derivation of the Geneva method has been presented extensively in several publications, e.g. in refs. [4,9]. Here, we explicitly refrain from entering into the finer details of the method, and we only briefly recall the general formulae, highlighting some key features that are important for this process. We use N -jettiness [68] to resolve the QCD emissions that can be associated with each event produced by Geneva: T 0 as the zero-jet resolution parameter, and T 1 to separate between one or more emissions. The partonic event space is then divided into three regions: Φ 0 for events with no extra emissions, Φ 1 for one-jet events, and Φ 2 for the remaining events with two jets in the final state. These phase space regions are defined via two thresholds, T cut 0 and T cut 1 . The differential cross section for the production of events with no extra emissions is given by Here we use the primed counting for the resummation order as in e.g. ref. [53]. For the case of a single extra emission we have two contributions: that above T cut 0 dσ mc 2) and the nonsingular below T cut 0 , arising from non-projectable configurations, Similarly the case of two extra emissions also receives two contributions, above and below T cut 1 , respectively. In the formulae above, B n , V n and W n are the 0-, 1-and 2-loop matrix elements for n QCD partons in the final state (including parton densities); analogously, we denote by N k LO n a quantity with n additional partons in the final state computed at N k LO accuracy. Since it is necessary to evaluate the resummed and resummed-expanded terms on phase space points resulting from a projection from a higher to a lower multiplicity, we introduce a shorthand for such projected phase space points, Φ N . We use the abbreviation to indicate an integration over the portion of the Φ M phase space which can be reached from a Φ N point while keeping some observable O also fixed, with N < M . The Θ O (Φ M ) term additionally limits the integration to the phase space points belonging to the singular contribution for the given observable O. For example, when generating 1-body events we use where the 1 → 2 mapping has been constructed to preserve T 0 , i.e.
and Θ T (Φ 2 ) guarantees that the Φ 2 point is reached from a genuine QCD splitting of the Φ 1 point. The use of a T 0 -preserving mapping is necessary to ensure that the pointwise singular T 0 dependence is alike among all terms in eqs. (2.2) and (2.4) and that the cancellation of said singular terms is guaranteed on an event-by-event basis. The non-projectable regions of Φ 1 and Φ 2 , on the other hand, are assigned to the cross sections in eqs. (2.3) and (2.5). These events are entirely nonsingular in nature. We denote the constraints due to the choice of map by Θ map , using the FKS map [69] for the Φ 1 → Φ 0 projection and, as mentioned above, a T 0 -preserving map for the Φ 2 → Φ 1 projection. Their complements are denoted by Θ map .
The term V C 1 denotes the contributions of soft and collinear origins in a standard NLO local subtraction, with C 2 a singular approximant of B 2 ; in practice we use the subtraction counterterms which we integrate over the radiation variables dΦ 2 /dΦ C 1 using the singular limit C of the phase space mapping.
In the formulae involving one or two extra emissions, U 1 is a next-to-leading-logarithmic (NLL) Sudakov factor which resums large logarithms of T 1 , and U ′ 1 its derivative with respect to T 1 ; the O(α s ) expansions of these quantities are denoted by U We extend the differential dependence of the resummed terms from the N -jet to the (N +1)-jet phase space using a normalised splitting probability P(Φ N +1 ) which satisfies The two extra variables are chosen to be an energy ratio z and an azimuthal angle ϕ. The functional forms of the P(Φ N +1 ) are in principle only constrained by eq. (2.10). However, in order to correctly model the soft-collinear limit behaviour, we find it useful to write them in terms of the Altarelli-Parisi splitting kernels, weighted by parton distribution functions (PDFs).
In previous implementations of the Geneva method, the splitting functions P(Φ N +1 ) were computed using a "hit-or-miss" method based on precomputed upper bounds, which did not require knowledge of an analytic expression for the integration limits of z and ϕ (see section II.B.4 of ref. [4] for the definition of the splitting function and Appendices C and D of ref. [70] for the computation of the upper bounds in a similar situation). At the same time, however, this introduced some numerical instabilities. In this work, we improve on this situation by including the exact integration limits and evaluate the splitting functions directly for each phase space point, as detailed in sec. 3.1.

Higgs boson production via gluon fusion
We consider the production of a stable Higgs boson via the gluon fusion channel in protonproton scattering, pp → H + X, where X denotes any additional hadronic radiation in the final state. At leading order (LO) in the strong coupling this results in a single contribution gg → H at partonic level [71], while at next-to-leading order (NLO) (anti)quark-initiated channels also start to contribute [72][73][74].
For a stable Higgs boson it is phenomenologically reasonable to work in the HTL effective field theory (EFT), in which the contributions from the top-quark loops coupling the Higgs boson to gluons have been integrated out. This EFT supplements the Standard Model (SM) vertices with additional, effective couplings between gluons and Higgs bosons. Introducing these effective vertices has the advantage of reducing the complexity of the matrix element computations. The cross section dependence on the top-quark mass m t can be partially restored by rescaling the HTL result by a factor equal to the ratio between the LO m t -exact result and that obtained in pure EFT. This is later referred to as rescaled EFT (rEFT), and reproduces the exact m t dependence of the LO cross section by construction. It is known to be a good approximation, for inclusive quantities, at least up to NNLO [64]. The resulting approximation can instead be problematic for differential distributions, for instance the transverse momentum of the Higgs boson when the accompanying radiation resolves the top-quark loop, i.e. when its transverse momentum is larger than m t . For the case of a finite top-quark mass, the NNLO corrections have been recently calculated for the inclusive cross section [63,64], and those at NLO for Higgs boson production in association with up to two hard jets [75,76]. At this level of precision, however, one also needs to take into account the interference between contributions including both massive top and bottom quarks, which is known at NLO for the Higgs plus jet case [77,78]. Since the problem of including the quark mass effects for precise phenomenological studies is largely independent of the matching of fixed-order and resummed calculations to parton showers in the Geneva method, which is the topic of the present study, we leave the investigation of these effects to future work.
In this work, the Higgs boson is always produced on shell with a mass m H = 125.09 GeV. When comparing with data in the fiducial regions of the ATLAS or CMS experiments, we will consider Higgs boson decays. In this case we work in the narrow-width approximation, which for the Higgs boson is particularly accurate since Γ H /m H ∼ O(10 −5 ). The Higgs decay products can always be added a posteriori due to the scalar nature of the boson, which implies that they are isotropically distributed without spin correlations with the initial state. For the rest of this work we will consider a collider energy of √ S = 13 TeV and assume the following values for the SM parameters affecting our calculations: (2.11) For the matrix elements in the HTL approximation we use the heftpphj and heftpphjj libraries of OpenLoops2 [79][80][81], which we then rescale by the rEFT factor r EFT = 1.06545.

T 0 resummation
The formulae presented in sec. 2.1 require the evaluation of the resummed spectrum and cumulant in the resolution variable T 0 up to at least NNLL ′ accuracy. Although the Geneva method does not depend on any particular resummation formalism, in practice we often find it convenient to exploit results derived via SCET [82][83][84][85][86][87][88]. Within this framework, a factorisation theorem for the zero-jettiness was first derived in refs. [89,90] for colour singlet production. In the case of the gluon-fusion channel for Higgs production it reads where H gg→H , S gg , and B g are the hard, soft and beam functions, respectively. The process-specific hard function H gg→H (Q 2 , µ) is defined as the square of the Wilson coefficient that results from matching the QCD Hamiltonian to the SCET operators, and encodes information about the Born and virtual squared matrix elements. It depends only on the Higgs boson virtuality Q 2 . In this section and whenever we consider Higgs boson production specifically, we set Q = m H ; elsewhere, we consider Q to be a generic hard scale.
Given that we work in the HTL approximation, we perform a two-step matching procedure: we first integrate out the hard degrees of freedom above the top-quark mass, and subsequently match the resulting EFT onto SCET. The final hard function then arises from the product of two Wilson coefficients, the first from the HTL approximation and the second from the matching to SCET; we evaluate both at the same scale µ. In principle, within this approach one could resum ln(m t /m H ) contributions by renormalisation group equation (RGE) evolution. However, given the values of the top quark and Higgs boson masses, these logarithms are never large and, consequently, we include them only at fixed order in the hard function. Alternatively, if one wants to include the full top-quark mass effects, a single-step matching can be performed as in e.g. ref. [53] at NNLL. Extending this to NNLL ′ accuracy requires the three-loop hard function with the exact top-quark mass dependence [63,64].
The beam functions B g (t, x, µ) are the inclusive gluon beam functions [89], which depend on the transverse virtualities t a,b of the initial-state partons that participate in the hard interaction and on their momentum fractions x a,b . While they are nonperturbative objects, for t ≫ Λ QCD they admit an operator product expansion (OPE), where the I ij (t, z, µ) are matching coefficients that describe the collinear virtual and real initial-state radiation (ISR) and the f j (ξ, µ) are the usual PDFs. For later use, we denote the Mellin convolution via the symbol ⊗ x . Finally, S gg (k, µ) is the gluon hemisphere soft function for beam thrust. Like the beam functions, S gg (k, µ) is a nonperturbative object and for k ≫ Λ QCD it also satisfies an OPE, where the LO matching coefficient is calculable in perturbation theory. Its perturbative component depends only on the colour representation of the hard partons, and therefore the gluon case can be derived from that of the quark channel via Casimir scaling. In our calculation we neglect the nonperturbative part of the soft function. We then rely on the hadronisation model of the parton shower to provide the missing contribution.
The functions in eq. (2.12) are all evaluated at a common scale µ and satisfy RGEs. The scale dependence in each of these functions involves potentially large logarithms of ratios of disparate scales, which may impact their perturbative convergence. In order to reduce the effect of these large logarithms, we evaluate each function at its characteristic (canonical) scale, i.e. µ S = T 0 , µ H = m H , and µ B = √ µ S µ H . Since the cross section needs to be evaluated at a common scale µ, we use the RGEs to evolve each function to µ. In doing so, we resum said logarithms at all orders in perturbation theory. The resummed formula for the T 0 spectrum is then given by (see e.g. ref. [9] for more details) where we denote the standard convolutions between the different functions and the RGE evolution factors via the ⊗ symbol. In order to achieve NNLL ′ accuracy in the T 0 resummation, each of the hard, soft and beam function boundary terms must be known at 2-loop order. For the beam function they were calculated at 2-loops in ref. [91], and in fact they are known up to 3-loop order [92]. Our implementation of the gluon beam function relies on an interface to scetlib [43,93,94], a library which provides ingredients for resummed calculations in SCET. The soft function has been known at 2-loops for some time [95,96], and recent work has aimed to push this calculation to the 3-loop order [97][98][99]. The hard function has appeared several times in the literature, see e.g. refs. [53,100], and is known analytically with full top-quark mass dependence at NNLO [63]. In addition, the anomalous dimensions and the beta function which enter the evolution factors and the fixed-order expansion of eq. (2.14) must be known at 2-loop (noncusp [53]) and 3-loop (cusp [101][102][103], β(α s ) [104,105]) order. By including them at one order higher [53,106,107], one can achieve resummation at N 3 LL.
The resummation of T 0 for the case of Higgs boson production via gluon fusion has already been studied in ref. [53] up to NNLL accuracy. In the present work, we extend this calculation to NNLL ′ and N 3 LL. For the determination of the canonical scales we employ the T 0 -dependent profile functions described e.g. in sec. 3 of ref. [9] with {x 0 , x 1 , x 2 , x 3 } = {1.5 GeV/m H , 0.2, 0.35, 0.5}. The use of such T 0 dependent scales is known to cause a difference between the integrated spectrum and the cumulant, which is formally of higher order. This is a result of the noncommutativity of the scale setting and the integration steps. In previous Geneva implementations, this problem has been alleviated by explicitly adding higher order terms to restore the cumulant cross section (see eq. (45) of ref. [9]). This can be done either by using a 'brute-force' approach, in which the integrated spectrum is simply replaced by the cumulant, or by smoothly transitioning from one to the other as a function of T 0 . In all Geneva implementations thus far we have followed the latter approach, which has the advantage of preserving the T 0 spectrum in its peak region.
In the case of gg → H production at 13 TeV, the difference between the integrated spectrum and the cumulant amounts to ∼ 18% of the total cross section. Given the size of these corrections, we found the previously adopted solution to be insufficient to completely solve the mismatch. In particular, our smooth fix modifies the T 0 spectrum in the region between ∼ 10 and ∼ 30 GeV by too large an amount, moving the central value of the first outside the uncertainty bands of the second. We therefore revert to the brute-force approach, and only require the preservation of the resummed cumulant cross section by fixing κ(T 0 ) = 1 (see eq. (45) of ref. [9]) such that the spectrum is exactly equal to the derivative of the cumulant.

Novel features of the GENEVA method
In this section we discuss the new improvements that have been incorporated in the Geneva method. Here we focus on their impact on the gg → H process, however we note that they can be straightforwardly generalised to several other processes (and indeed have already been tested for Drell-Yan, double Higgs [15], and tt production [16]).

Improved treatment of splitting functions
The N -jettiness spectra computed through resummation techniques cannot be directly used for generating events with N + 1 final-state partons, since they do not carry a dependence on the full Φ N +1 configurations, but only on T N and the projected Φ N configurations with N final-state partons. For this reason, a splitting function P(Φ N +1 ) was introduced in ref. [8] in order to make the resummed calculation fully differential in the higher order phase space.
In general, the N → N + 1 splitting function If the function g(Φ N , T N ) is the T N spectrum, then multiplying it by the P functions makes it differential over the dΦ N +1 phase space without affecting the distributions of observables that only depend on Φ N and T N . In order to provide an explicit expression for P, we write the phase space of the Φ N +1 configurations with a valid Φ N projection as the product of dΦ N , dT N and the phase space parametrised by two additional radiation variables z and ϕ. In this way the integral over the projectable Φ N +1 configurations at fixed Φ N and T N can be expressed as where the index k runs over the N + 2 possible emitter partons (mothers) of the Φ N configurations. For each mother k and its associated mapping, we assume that the Jacobian does not depend on ϕ. This is true for all the mappings considered in this paper. The integral over the Φ N +1 configurations summed over the n real partonic subprocesses with N + 1 final-state partons for a generic function g β (Φ N +1 ) can now be written as where n Born is the number of subprocesses with N final-state partons, and n split k the number of possible QCD splittings k → i+j, with i the emitted parton and j the sister. The function g k→i+j α (Φ N , T N , z, ϕ) on the right-hand side is equal to g β (Φ N +1 ) expressed in terms of the underlying Born process index α and the splitting indices k and j. For ease of notation, the full dependence of the z and ϕ integration limits on the phase space variables is not shown. The unprojectable Φ N +1 configurations are those for which either the two closest partons do not represent a valid QCD splitting, the Φ N configuration obtained from the projection is not kinematically allowed, or the flavour configuration of the Φ N is invalid.
In order to fulfil the condition presented in eq. (3.1), we choose splitting functions P(Φ N +1 ) that depend on the mother and sister indices and vanish in the unprojectable Φ N +1 configurations: The P kj must then satisfy the equation for all values of Φ N and T N . Without loss of generality, in the projectable Φ N +1 configurations we can express them as where f kj is a generic function that we specify later. If we choose it to be independent of ϕ, the above expression simplifies to . In order to perform the integral in the denominator of eq. (3.8), we compute the integration limits on z and ϕ and the Jacobian J k both for the 0 → 1 and 1 → 2 splitting mappings for each Φ N +1 configuration. In the previous Geneva implementation of the splitting functions the computation of the integration limits was avoided by precomputing their upper bounds and then using a "hit-or-miss" integration method. We highlight that, whenever the constraints on z and ϕ are in the form of an inequality involving both the variables, we only compute an overestimate of the true integration limits on z analytically. We then determine the true limits numerically by imposing the condition ∆ϕ k (Φ N , T N , z) > 0.

Infrared limits
For this section we introduce the acronyms ISRA (initial-state radiation A), ISRB (initialstate radiation B), and FSR (final-state radiation) to indicate the N + 2 possible mothers we have to deal with: the parton from the first (A) and second (B) beam, respectively, and the final-state partons. We furthermore label ISRA and ISRB collectively as ISR.
The exact form of the function f kj in eq. (3.7) can significantly affect the efficiency of the Monte Carlo event generator in the region of small T N > T cut N . In this region, the logarithmically enhanced terms appearing in the fixed-order calculation have to cancel those coming from the resummed-expanded contributions. For this reason the main criterion we follow in the choice of f kj is to achieve a good approximation of the behaviour of the associated matrix element in the infrared limit when T N → 0.
For simplicity, in practical applications we choose not to include the azimuthal dependence in the form of the f kj functions, using eq. (3.8). We define where a and b are the initial-state partons, α s (µ R ) is the strong coupling evaluated at the renormalisation scale µ R , and f H i (x i , µ F ) is the PDF of the parton i in the hadron H evaluated at longitudinal momentum fraction x i and factorisation scale µ F . The renormalisation and factorisation scales are fixed to µ R = µ F = Q, where Q is the virtuality of the colour singlet system. TheP kj are the unregulated Altarelli-Parisi splitting functionŝ We highlight that for the 0 → 1 splitting, connecting events with no extra partons to events with one extra parton, the PDFs are evaluated at the exact momentum fractions x a (z) and x b (z) of the real emission phase space Φ 1 rather than their infrared limits. This has proven to be necessary to obtain an accurate description also in the tail of the colour singlet transverse momentum distribution. We note that in this case we also reproduce the correct soft limit, as shown in sec. 3.1.2. For the 1 → 2 splitting the true x a and x b also depend on ϕ. In this case they are approximated by dropping this additional dependence, which still represents an improvement with respect to the strict collinear limit.

Soft limit of the 0 → 1 splitting
In the following we show that the expression of f kj introduced in eq. (3.9) correctly reproduces both the singular soft and collinear limits at O(α s ) in the 0 → 1 splitting.
In the case of colour singlet production in hadron-hadron collisions, let us consider the k → i + j splitting connecting the Born matrix element B 0 and the real matrix element B 1 (in both cases excluding parton densities). This can be expressed in terms of the FKS variables ξ = 2 E/ √ s = 1 − Q 2 /s and y = cos θ. Here, s is the partonic centre-of-mass energy squared, and E and θ are the energy of the emitted parton and the angle between the emitted and the right-moving incoming parton in the partonic centre-of-mass frame, respectively.
In the soft limit of the emitted particle i, we have where C k = C F for the quark-initiated processes and C k = C A for the gluon-initiated. In the azimuthally averaged collinear limit between particles i and j, we have where y → 1 and y → −1 represent the collinear limits with respect to incoming parton a and b respectively. If the colour singlet production process is quark-initiated or has only scalar particles in the final state, the above expressions also hold prior to averaging over the azimuthal angle.
We consider a configuration with one final-state parton with momentum p, where Herep is obtained by longitudinally boosting p from the laboratory frame to the frame where the colour singlet has zero rapidity, andp ± =p 0 ∓p 3 . We have chosen z such that in the collinear limit it reduces to the energy fraction of the emitter with respect to the sister, while providing the correct scaling also for the single soft limit.
In order to show that the singular limits in eq. (3.9) reproduce the above results, we rewrite T 0 and z in terms of the FKS variables ξ and y, and then compare the ensuing expression to eqs. (3.11) and (3.12). They read (3.14) Therefore in the infrared singular limit one obtains

16)
z → 1 − ξ in the collinear limit. (3.17) Multiplying the NLL singular T 0 spectrum expanded at O(α s ) by the splitting functions, in the infrared limit we find up to power corrections. By using the above expressions for T 0 and z, it can be shown that this reproduces both the soft and collinear limits given in eqs. (3.11) and (3.12). We remark that with the choice of z given in eq. (3.13) the soft limit can be entirely captured by using the Altarelli-Parisi splitting collinear kernels. The validity of eq. (3.18) can be understood to be a consequence of the fact that the noncusp soft anomalous dimension is zero at one loop order, resulting in the lack of a single logarithmic contribution to the T 0 spectrum coming from the soft function.

Numerical validation
In this section we present the effects of the improved splitting function P impr implementation described above in the case of Higgs boson production via gluon fusion, setting Q = m H ; we focus on the p H T and the T 0 spectrum. We compare the results of a fixedorder calculation with those obtained by truncating the resummation formula in eq. (2.14) multiplied by the splitting function to the same order. We do so for the results at LO 1 compared with the NLL ′ resummed-expanded in fig. 1, and for those at NLO 1 compared with the NNLL ′ resummed-expanded in fig. 2. We also show the nonsingular contribution, defined as the difference between these fixed-order and resummed-expanded pieces. In all plots we also show the results obtained with the original implementation P orig of the splitting function in eq. (3.8), which was based on a hit-or-miss method using upper bounds tabulated on a grid.
We begin the discussion with the results for the T 0 distribution. As expected, the improved implementation gives identical results to the original, both at LO 1 and NLO 1 . This is a consequence of the fact that T 0 is preserved by the splitting, by construction. We observe that at extremely low values of T 0 the presence of technical cuts in the fixed-order calculation affects the convergence to the singular predictions in both approaches. When instead considering the LO 1 results for the p H T distribution, we notice how the improved implementation of the splitting functions correctly captures the logarithmic behaviour of the matrix element at fixed order. This can be seen by the fact that the improved nonsingular distribution converges to zero, contrary to the original case which converges to a finite value. Similarly, an improvement is also visible for the NLO 1 case. Here, however, the new splitting function P impr is not able to exactly reproduce the complete logarith- mic behaviour of the NLO 1 result, as it appears to miss a single logarithmic contribution ∼ 1/p H T . This is implied by the fact that the improved nonsingular contribution converges to a nonzero constant at low values of p H T . This must however be compared with the original approach, P orig , where the divergent behaviour of the nonsingular plot suggests that that implementation also fails to capture the logarithmic structure up to ∼ ln 2 (p H T )/p H T . We examine the effects of the P impr implementation on the Drell-Yan process in App. A, where we compare different Geneva results with the ATLAS experimental data.

Independent scale variations
In traditional implementations of fixed-order QCD calculations, a differentiation is made between the factorisation scale µ F and the renormalisation scale µ R . The former is associated with the scale of collinear factorisation, while the latter is introduced in dimensional regularisation in order to render the strong coupling dimensionless.
To date, implementations of Geneva have assumed these scales to be equal. Doing so facilitated the matching to the resummed calculation, where a sole "nonsingular" scale µ NS appears as the endpoint of the RGE running, typically taken to be a hard scale Q of the problem. The two scales were then varied in a correlated fashion ("diagonal" in the {µ R , µ F } space) when probing the higher order uncertainties. This approach, however, can hinder a complete and thorough uncertainty estimation as it neglects those variations which are off-diagonal, i.e. where µ R and µ F are varied independently. In this section we provide an improved and robust uncertainty estimation within the Geneva framework by exposing the µ F dependence of the singular cross section that eventually allows for off-diagonal scale variations, and discuss the choice of µ F in the infrared region.

Exposing the µ F dependence of the singular cross section
The collinear beam functions B i entering the T 0 factorisation in eq. (2.12) satisfy the OPE in eq. (2.13). In resummed predictions, they are evaluated at a scale µ = µ B where µ B minimises the singular logarithmic structure of B i , whereas at fixed order µ = µ R = µ F = Q, where for example Q = m H for on-shell Higgs boson production.
In order to expose the µ F dependence of the beam functions, we rewrite eq. (2.13) as where we dropped the power corrections. Here we evolve the PDFs from µ F to µ using the evolution kernel U ij (x, µ, µ F ) that results from the solution of the DGLAP equations, and we follow the conventions of ref. [93] for the perturbative expansion of the splitting kernels P ij (x, µ). 1 Although the µ F dependence in eq. (3.19) cancels exactly between the PDFs and the evolution kernel, as soon asÎ ij is truncated at a given order, a residual µ F dependence appears in the beam function, In order to manifest this dependence explicitly, we note that the matching coefficients and the evolution kernel in eq. (3.19) admit perturbative expansions in the strong coupling constant,Î We first obtain closed-form expressions for the U ij , which can be achieved either by directly solving eq. (3.20) or by using the solutions of the RGE satisfied by I ij (see eq. (2.17) of ref. [93]). We have where U −1 ij denotes the inverse of U ij . Here and in the following, repeated flavour indices are implicitly summed over. We note that eq. (3.24) is exactly the same as eq. (2.17) of ref. [93] if we set γ B = γ ν = 0. It is therefore straightforward to use its solution, which up to O(α 2 s ) reads where we have abbreviated L = ln(µ/µ F ). Solving the closure equation of the evolution kernels where the expressions for the (µ F -independent) matching coefficients I ij (t, x, µ), the anomalous dimensions γ i B 0 and Γ i 0 , and the plus distributions L n can be found in ref. [93].

Choice of the factorisation scale
The choice of the factorisation scale µ F is in principle subject to different requirements in the fixed-order and in the resummation region. In order to minimise the size of the logarithms L = ln(µ/µ F ) in eqs. (3.30) and (3.31), in the resummation region we demand that µ F ∼ µ B , i.e. the beam scale. On the other hand, in the fixed-order region, a natural scale setting is µ F ∼ µ R ∼ µ NS so that the fixed-order perturbative convergence is not jeopardised. However, given that the beam function profile scale µ B (T 0 ) flows to µ NS in the fixed-order region, choosing µ F = µ B satisfies both conditions. When considering the scale variations for the estimation of the theoretical uncertainty, the definition of the profile scales is extended to where the central predictions are obtained setting κ R = κ F = 1, α = α f = β = 0. The function f run is defined in ref. [9] and the function f vary in ref. [108]. We note that a different parameter is used in the exponent of the function f vary for µ B and µ F in eqs. (3.34) and (3.35) above, while we use the same β parameter for both. This is justified because the β-variations are introduced in order to disentangle the variations of µ B and µ S ; no ratios of µ F /µ S appear in the singular cross section, therefore there is no need for an independent β parameter in µ F . We also extend the fixed-order uncertainty ∆ FO to include the off-diagonal µ F and µ R = κ R µ NS scale variations, resulting from the envelope of a 7-point scale variation

Treatment of timelike logarithms
Radiative corrections in colour singlet production processes such as Drell-Yan or gluoninitiated Higgs production contain Sudakov logarithms of the form α n s ln m (−q 2 /µ 2 R ) , m ≤ 2n, where q µ is the momentum of the colour singlet system. They primarily arise in the calculation of the corresponding form factor, and their coefficients are linked to the structure of its infrared singularities [109]. For such processes q µ is a timelike vector, i.e. q 2 > 0, and the scale choice µ 2 R = q 2 results in the Sudakov logarithms developing an imaginary part, since ln m (−1) = (±iπ) m . The presence of such 'timelike' logarithms at each order might negatively affect the perturbative convergence of the cross section, where they result in additional terms proportional to powers of π 2 . The severity of this effect is process specific; for Drell-Yan or gluon-initiated Higgs production it has been explicitly studied for both exclusive [45,46,53,55,90,108,[110][111][112] and inclusive [19,43] observables. A way to mitigate the impact of said contributions is to choose µ R such that timelike logarithms are eliminated, i.e. to evaluate the form factor at the complex scale µ R = −i|q| = −i Q [113][114][115][116].
In factorised singular cross sections, the square of the form factor naturally appears as the hard function H(Q 2 , µ H ) of the process. Following the above discussion, the hard function can be evaluated at the complex scale µ H = −i Q. This choice implies a nontrivial renormalisation group evolution between µ H and the real scale Q, in the form of a rotation in the complex µ plane. This procedure is also referred to as 'timelike resummation'. It has been applied to a multitude of exclusive and inclusive processes, including the case of Higgs boson production, see e.g. refs. [18,19,43,53,90,117].
In addition, it has been shown that not only the singular, but also the nonsingular [53,55] and the total [19] cross section can benefit from this prescription, in terms of an improved perturbative convergence and a reduced ensuing scale uncertainty. This can be better understood when one considers that by integrating the timelike resummed exclusive cross section, one obtains the corresponding resummed inclusive prediction. For the nonsingular contribution, a factorisation formula remains in general unknown, meaning that it cannot be directly evaluated at the complex scale µ H -it still, however, contains the form factor plagued by timelike logarithms. In this case, the procedure for the treatment of these logarithms involves re-expanding the nonsingular contribution, extracting from it the hard function evaluated at Q, and replacing it with that evaluated at the scale µ H , as detailed in ref. [19].
In our implementation we perform these steps at the same order as the T 0 resummation, both in the singular and the nonsingular terms. This implies that the improved perturbative convergence following the choice of a complex-valued scale µ H = −i Q does not only apply to the singular T 0 spectrum but also to the inclusive predictions.
In order to study the uncertainty associated with this choice of scale, we follow the prescription introduced in ref. [19], designed to probe the structure of the timelike logarithms. The uncertainty ∆ φ is estimated by the envelope of the phase variations while the central value predictions correspond to φ = π/2. Since there is no dynamical parameter governing the choice of the scale µ H , the timelike resummation is performed throughout the T 0 spectrum, i.e. even when T 0 resummation is off. We therefore consider the uncertainty resulting from variations in eq. (3.38) as an independent source and add it to the other uncertainties in quadrature. Thus, for inclusive predictions we have whereas for exclusive predictions we use In the Geneva implementation of the gg → H process we use the hardfunc module from scetlib [94] for the hard function evaluation and evolution in the complex plane. Since for this process we set Q = m H , we pick With this choice, we observe a difference in the total cross section result with respect to the µ H = m H case that can be substantial despite being formally of higher order. The effects of the complex choice of scale µ H on differential observables are illustrated in fig. 3, where we compare predictions at NNLO+NNLL ′ for the T 0 and y H distributions with µ H = m H and µ H = −i m H . In this and the following figures, the theoretical uncertainty is shown as a shaded band, while the Monte Carlo integration errors are shown as thin vertical bars. For the Higgs boson rapidity distribution, we observe an increase of around 10% that is almost independent of y H , and a reduction in the uncertainty band as expected. The T 0 spectrum shows a larger effect, especially in the tail of the distribution, where our prediction is entirely driven by the fixed-order result. Nonetheless, we observe a reduction in the uncertainty band particularly in the peak and transition regions of the spectrum, between 5 and 45 GeV.

Validation of the gg → H process
In this section we validate our predictions. We first compare our partonic NNLO results with two independent calculations, and then discuss the interface to the Pythia8 shower.

Partonic results at NNLO
Here we validate the NNLO accuracy of the total cross section obtained with Geneva and that of the only differential inclusive quantity available, the Higgs boson rapidity. We compare the total cross section with the independent calculations implemented in ggHiggs [62, [118][119][120][121] and Matrix [122], and the rapidity distribution with Matrix only. The Matrix predictions are based on the q T -subtraction approach and are extrapolated towards the zero q T -cut value. We set the input parameters of our calculations as described in sec. 2.2, and we choose the central factorisation and renormalisation scales equal to each other and to the Higgs boson mass, µ R = µ F = m H . We set our resolution cutoffs to T cut 0 = T cut 1 = 1 GeV. We employ the PDF set PDF4LHC15 nnlo 100 from LHAPDF [123], and take the value of α s (m Z ) from the same set, so that α s (m H ) = 0.11263.
In table 1 we report the values of the inclusive gg → H cross section and the associated 7-point scale variations calculated at NNLO and rescaled with the rEFT factor using Geneva, ggHiggs, and Matrix. 2 We observe excellent agreement between the three predictions; by choosing T cut 0 = 1 GeV, the neglected power-suppressed terms in Geneva are at the permille level and amount to an acceptable ∼ 0.02 pb error for the central value.
In fig. 4 we compare the Higgs rapidity spectrum obtained with Geneva with the NNLO result provided by Matrix, including the 7-point scale variations. We observe very good agreement both in the central values and in the envelope of the scale variations, up to large values of |y H |. The symmetry of the pp collider allows us to show only the absolute value of y H , and thus further reduce the Monte Carlo uncertainty.

Interface with PYTHIA8
In this section we briefly recap the main features of the interface used in Geneva to match the partonic results to the Pythia8 [124] parton shower. As this is not the main focus of this work, however, we refer the interested reader to ref. [4] for a detailed discussion and ref. [15] for additional details on the accuracy of the matched calculation. Given that so far we have constructed partonic results with NNLL ′ accuracy in the resolution variable T 0 , we wish to preserve this resummed accuracy after the parton shower as far as is possible. At the same time, for all other observables we need to guarantee that the accuracy of the parton shower is preserved. This is a nontrivial condition: since the ordering variable of the Pythia8 parton shower is the relative transverse momentum while the resolution variable we use is the N -jettiness, the shower can in principle produce emissions which double-count regions of the phase space.
To avoid this issue, we perform the matching employing the following prescription. We set the starting scale of the parton shower by taking the maximum relative k ⊥ determined by the lower scale of the resummation. The latter is defined on an event-by-event basis and corresponds to either T c N ≡ T cut 0 , T cut 1 or T 1 (Φ 2 ), depending on whether the relative partonic configuration has N = 0, 1 or 2 jets, respectively. We then let the shower run down to the internal minimum p ⊥ , which produces a certain number of emissions k. Lastly, we check that the resulting event fulfils the condition which ensures that both accuracies are correctly preserved. For unshowered events with one jet in the final state, we perform the first shower emission directly within Geneva, by implementing eqs. (48) and (49) of ref. [9]. Showered events will therefore almost exclusively originate from events with either zero or two final state partons.
In fig. 5 we show the effect of the Pythia8 shower on the p H T and y H partonic distributions. For the results presented in this section we use the default Pythia8 parameters for the shower and the hadronisation model. The rapidity distribution, being an inclusive observable, is exactly preserved by the shower, as expected. The Higgs transverse momentum is an exclusive observable, and the shower can therefore have a significant impact on its shape: in this case we see an effect of ∼ 15% in the p H T < 15 GeV bin, and smaller effects ≲ 5% in the rest of the spectrum, especially in the tail of the distribution. After hadronisation, we find that most of these discrepancies are reduced. The parton shower and hadronisation effects on the T 0 distribution are displayed in fig. 6. As mentioned above, our matching procedure to the Pythia8 shower is designed with the aim that the T 0 logarithmic accuracy is not spoiled. We explicitly check this in the left panel, where we compare the T 0 distribution at NNLL ′ before and after the parton shower matching with the partonic prediction at N 3 LL. Note that for this process, which is gluon-initiated, one expects that the parton shower effects are larger than for quark-initiated processes, e.g. because of the larger Casimir factors. Nonetheless, in the peak region T 0 < 25 GeV, we find that the showered distribution lies in between the central NNLL ′ and N 3 LL curves, and within the overlap of the two uncertainty bands. We therefore conclude that the quantitative effects of the shower are on par with (or smaller than) the effects of the next logarithmic order in the resummation.
In principle it is possible to directly interface the N 3 LL Geneva results to the parton shower. In this work, however, we refrain from doing so -in particular when comparing to data -because, due to the lack of the N 4 LL prediction for the T 0 spectrum, we cannot verify that the large distortions induced by the shower are compatible with the next logarithmic correction.
The hadronisation effects on the T 0 distribution are displayed in the right panel of fig. 6. As expected for this observable, we observe O(1) effects in the peak region, which decrease for larger values of T 0 . In the region around T 0 ≈ m H /2, which corresponds to the point at which the T 0 resummation is switched off, we find a more pronounced discrepancy between the Geneva partonic and showered results. We have verified that this is an artefact related to our choice of setting the T 0 spectrum equal to the derivative of the cumulant as explained at the end of sec. 2.3.

Comparison with LHC data
We compare the predictions obtained with Geneva with the latest experimental results for the Higgs boson inclusive and differential cross sections in the H → γγ decay channel. The results are provided both by the ATLAS [66] and CMS [67] experiments, and are obtained from the LHC data at a centre-of-mass energy of 13 TeV using 139 fb −1 and 137 fb −1 of proton-proton collision data, respectively.
In the ATLAS measurement [66], the fiducial phase space is identified by requiring the existence of two isolated photons with p γ T > 25 GeV in the final state. Photons are considered isolated if the transverse energy of charged particles with p T > 1 GeV within a cone of radius R iso = 0.2 around the photon direction does not exceed 5% of the photon's transverse momentum. The two isolated photons must additionally have transverse momenta larger than 35% and 25% of the diphoton invariant mass, for the leading and subleading photons respectively. The invariant mass of the diphoton system must be in the range 105 GeV < m γγ < 160 GeV. Moreover, photons are required to have pseudorapidity |η γ | < 1.37 or 1.52 < |η γ | < 2.37. For this measurement, jets are defined using the anti-k T algorithm with radius R = 0.4, and must have p j T > 30 GeV and |y j | < 4.4. Jets must also be separated from photons with p γ T > 15 GeV by a distance ∆R γj > 0.4.
Similarly, in the CMS measurement [67], the fiducial region is defined by having two isolated photons in the final state. In this case, photons are isolated if the transverse energy of all particles inside a cone of radius R iso = 0.3 is less than 10 GeV. The transverse momenta of the leading (subleading) isolated photon must satisfy p γ T > 35 (25) GeV, and amount to at least 1/3 (1/4) of the reconstructed Higgs invariant mass. In turn, the Higgs invariant mass must lie between 100 and 180 GeV. Photons must also satisfy |η γ | < 2.5. Also in this case, jets are constructed using the anti-k T algorithm with R = 0.4, and are required to have p j T > 30 GeV. Jets with |η j | < 2.5 are used for observables with one extra jet or to count the number of jets, while a looser cut |η j | < 4.7 is applied for observables requiring at least two jets in the final state.
Due to the lack of availability of these analyses in the Rivet [125] framework, we have implemented the ATLAS and CMS analyses within the Geneva code. The H → γγ decay is inserted by the Pythia8 particle decays handler on top of the events produced by Geneva. Its kinematics are treated at leading order in QCD, and we set the branching ratio to BR(H → γγ) = 2.27 × 10 −3 , i.e. the value reported in ref.
[126] and calculated with HDECAY [127]. The Geneva prediction for the gluon-fusion production channel is obtained at NNLO+NNLL ′ T 0 +NLL T 1 , and setting the scale of the hard function µ H = −i m H . We use matrix elements computed in the infinite-top-mass limit and rescaled in the rEFT scheme. We set T cut 0 = T cut 1 = 1 GeV. We use the PDF set PDF4LHC15 NNLO, and take the value of α s (m Z ) from there. The partonic prediction is matched to the Pythia8 QCD+QED shower, including multiparton interaction (MPI) contributions. We use the AZNLO tune [128] for the ATLAS comparison, and the CP5 tune [129] for CMS. Showered events are then hadronised using the default Pythia8 Lund string model [130,131]. In order to obtain a meaningful comparison with the experimental data, we include the contributions from other Higgs boson production modes (labelled overall as XH) by summing them to the Geneva results for the gluon-fusion channel alone. 3 For ATLAS these include vector-boson fusion (VBF), Higgsstrahlung (V H), and associated production with tt, bb, and t, all computed at NLO accuracy in QCD. For CMS these only include contributions from VBF, V H, and ttH. The outcome of the comparison with the experimental results is shown in figs. 7 and 8 for the ATLAS data, and in figs. 9 and 10 for the CMS data. For the ATLAS data, we show the p H T , N jets , |y H |, p j 1 T , p Hj T , and τ j 1 C distributions, as well as the p H T spectra in bins of τ j 1 C and in bins of |y H |. For the CMS data, we show the p H T , N jets , |y H |, p j 1 T , τ j 1 C , and |ϕ * η | distributions, as well as the p H T spectra in different jet multiplicity bins (N = 0 and N = 1). The definitions of τ j C and ϕ * η are given by where ϕ acop = π − |∆ϕ γγ | and sin θ * η = [cosh(∆η γγ /2)] −1 [132].
With the ATLAS fiducial cuts, we obtain a total fiducial cross section of 58.8 +1.5 −3.0 fb, to be compared to the experimental finding of 67 ± 6 fb. In the CMS fiducial region, we obtain a total cross section of 66.6 +1.6 −3.3 fb, which is compatible with the measurement of 73.4 +6.1 −5.9 fb. In both cases our predictions agree with the measured results within roughly one standard deviation. We note that our results are not rescaled to the total N 3 LO gluonfusion cross section, contrary to the theoretical predictions used in the ATLAS publication for their comparison.
Regarding the distributions, we find overall good agreement between the Geneva predictions and the measurements. For the ATLAS data we find slight deviations in the p H T peak and a more marked discrepancy in the tail of the distribution. The latter corresponds to the region where the HTL approximation is less accurate. We also find slight deviations in the |y H | spectrum. The deviations in both spectra are consistent with those obtained using other calculations, as shown in ref. [66]. Similarly, for the CMS data our results underestimate the bins corresponding to the p H T peak, again in a similar fashion to other predictions [67]. Large deviations are also found in the first bin of the p H T distribution with N jets = 1, once again in agreement with other theoretical predictions.

Conclusions
We have described a number of improvements to the Geneva method, which are particularly useful for all colour singlet production processes. Specifically, we detailed a new implementation of the splitting functions which serve to make the resummed calculation fully differential in higher multiplicity phase spaces. This results in an improved behaviour of the nonsingular cross section as a function of the colour singlet transverse momentum in the infrared limit. In addition, following earlier work [43] we have introduced a separation between the beam scale in our SCET-based resummed calculation µ B and the scale associated with collinear factorisation µ F which appears in the fixed-order calculation. This allowed us to achieve a more robust estimate of the theoretical uncertainties associated with our calculation in the fixed-order region, including uncorrelated variations of the renormalisation and factorisation scales. Finally, we have addressed the issue of large contributions from π 2 terms originating from timelike logarithms in 2 → 1 processes, by enabling the choice of a complex-valued hard scale µ H . We studied the associated resummation of said logarithms in our fully-differential calculation, and showed that, as previously noted in the literature, the perturbative convergence can thus be improved.
Throughout this work, we have used the gluon-initiated Higgs production process to study the effects of our improvements. We have constructed an NNLO+PS event generator for the process, including the resummation of the zero-jettiness variable up to NNLL ′ accuracy. We also studied the effects of the parton shower on the logarithmic accuracy achieved at partonic level by comparing the showered results to the N 3 LL partonic predictions. The availability of recent experimental results [66,67] for this process also allowed us to make a detailed comparison of our final, showered events with data. We stress, however, that the issues which we addressed in this work have a more general applicability, and we anticipate that the future implementations of processes in Geneva will make use of these developments.
In this study, we have consistently worked in a heavy-top limit in which the top-quark has been integrated out of the SM Lagrangian, resulting in an effective gluon-Higgs coupling. We reweighted the results with the exact LO top-quark mass dependence using the so-called rEFT approximation. Given the advancement in recent years towards including the exact quark mass dependence at NNLO [63,64], it would be desirable to incorporate this progress into a Geneva event generator at NNLO+PS. We leave this issue to future work.
The code used for this study will be included in a future public release of Geneva, and is available upon request to the authors together with the associated generated events.  Figure 11: Comparison of the ATLAS normalised p ℓℓ T and ϕ * η distributions [133] with the Geneva+Pythia8 results at 13 TeV. We show the distributions obtained with Geneva using T 0 as resolution variable and by using the original and improved splitting functions, and also the one obtained with Geneva+RadISH using p ℓℓ T as resolution variable.
predictions, we take the settings used in ref. [12]. We set T cut 0 = T cut 1 = p ℓℓ, cut T = 1 GeV, m Z = 91.1876 GeV, Γ Z = 2.4952 GeV, and α e (m Z ) = 7.55638 × 10 −3 . We use the NNPDF31 nnlo as 0118 PDF set, and take the value of α s from there. All the Geneva results are matched to the Pythia8 QCD shower including MPI contributions.
The results of this study are shown in fig. 11. With the improved splitting functions we get a p ℓℓ T distribution that is closer to data in the interval [10,100] GeV, although not as close as in the p ℓℓ T -resummed case. This is an effect of the better physical description of the splitting behaviour encoded in the expressions defined in eqs. (3.8) and (3.9). For the ϕ * η distribution, the results are affected by the parton shower to a greater extent, and the different Geneva implementations produce mixed performances, again with the p ℓℓ Tresummed case being the closest to data. We notice, though, that in the small ϕ * η limit using the improved splitting functions in the T 0 -resummed case gives a result that is closer to data than the one obtained with the original splitting functions. [23] ATLAS collaboration, Measurements of fiducial and differential cross sections for Higgs boson production in the diphoton decay channel at √ s = 8 TeV with ATLAS, JHEP 09 (2014) 112 [1407.4222].
[24] ATLAS collaboration, Fiducial and differential cross sections of Higgs boson production measured in the four-lepton decay channel in pp collisions at √ s=8 TeV with the ATLAS detector, Phys. Lett. B 738 (2014) 234 [1408.3226].
[25] ATLAS collaboration, Measurement of fiducial differential cross sections of gluon-fusion production of Higgs bosons decaying to WW * →eνµν with the ATLAS detector at [32] CMS collaboration, Measurement of inclusive and differential Higgs boson production cross sections in the diphoton decay channel in proton-proton collisions at √ s = 13 TeV, JHEP 01 (2019) 183 [1807.03825].
[33] CMS collaboration, Measurement and interpretation of differential cross sections for Higgs boson production at