Higgs boson decay into four leptons at NLOPS electroweak accuracy

In view of precision studies of the Higgs sector at the Run II of the LHC, the improvement of the accuracy of the theoretical prediction is becoming a pressing issue. In this framework, we detail a calculation of the full Next-to-Leading Order (NLO) electroweak corrections to Higgs boson decay into four charged leptons, by considering the gold-plated channel H ->Z(*) Z(*) ->2l 2l', l,l' = e, mu. We match the NLO corrections with a QED Parton Shower (PS), in order to simulate exclusive multiple photon emission and provide novel results at NLOPS electroweak accuracy. We compare our NLO predictions to those of the program Prophecy4f and present NLOPS phenomenological results relevant for Higgs physics studies, with particular attention to precision measurements of the Higgs boson mass, spin-parity assignment and tests of the Standard Model. Our calculation is implemented in a new code, Hto4l, which can be easily interfaced to any generator describing Higgs boson production. As an example, we provide illustrative results for Higgs production and decay in the process gg ->H ->4l using POWHEG with NLOPS accuracy in the production mode.


Introduction
With the announcement in 2012 of the discovery of a new particle in the search for the Standard Model (SM) Higgs boson by the ATLAS [1] and CMS [2] collaborations at the CERN LHC, particle physics entered a new era. The data collected at the centre-of-mass (c.m.) energies of 7 and 8 TeV have been analyzed by the two experiments in order to establish whether the newly discovered particle is actually the boson predicted in the SM as relic of the mechanism of electroweak symmetry breaking (EWSB) [3][4][5][6][7][8].
The mass of the observed particle has been precisely measured by studying the two cleanest decay channels given by the decays into a photon pair and into four charged leptons. The combination of the two channels H → γγ and H → 4 (4e, 4µ, 2e2µ), which have excellent mass resolution and where excesses with large significance are observed [9][10][11][12][13][14], presently provides a mass measurement of approximately 125 GeV for each experiment, with a relative uncertainty of better than 0.2% for the combined ATLAS-CMS measurement.
Concerning the main production mechanisms of the SM Higgs boson at hadron colliders, i.e. gluon-gluon fusion, vector boson fusion (VBF), associated production with a massive vector boson and associated production with top quarks, the studies performed at the LHC, based on the analysis of individual production signal strengths for various decay modes, have provided a clear observation of Higgs production through gluon fusion and an evidence for VBF production, with a significance above the 3σ level, and for associated V H(V = W, Z) production at about 3σ [12,15]. Various tests of the couplings of the new particle to bosons and fermions have been carried out both by ATLAS and CMS collaborations. In particular, the measured ratio of the couplings of the Higgs particle to W and Z bosons, which is an important probe of the EWSB mechanism as fixed by the custodial symmetry, is compatible with the SM expectation and, more generally, no significant deviation from the SM is observed from the coupling strength studies [9,12,15]. Noticeably, evidence for the direct coupling of the Higgs boson to down-type fermions has been reported through the study of the challenging decay modes Higgs into bottom quarks and τ leptons [16,17].
Last but not least, the spin and parity quantum numbers of the discovered particle have been assessed, by means of a systematic analysis of its production and decay processes. The data strongly favor the scalar nature J P = 0 + of the observed particle, while rejecting other non-standard hypotheses (J P = 0 − , 1 ± , 2 + ) or possibility of CP mixtures at high confidence [14,18,19].
All these measurements marked the start of a new era of precision Higgs physics and were accompanied by an impressive theoretical effort summarized in three CERN reports by the LHC Higgs Cross Section Working Group [20][21][22]. These studies, as well as the related theoretical work, are in continuous progress and will continue during the Run II of the LHC at higher energies and luminosity.
In this paper, we focus on the Higgs boson decay into four charged leptons, i.e. H → Z ( * ) Z ( * ) → 4 , in order to provide novel precision predictions of interest for future studies of the Higgs sector at the LHC. This decay channel plays a particularly relevant rôle, as it provides the cleanest experimental signature, given by a peak in the four lepton mass spectrum on top of a flat and small background. Actually, the H → 4 decay mode allows to derive a precise mass measurement in the different combinations of lepton final states, to assess the spin-parity quantum numbers using sensitive angular distributions and to perform precision tests of the SM at the level of differential cross sections [23]. In the off-shell region, the H → 4 data can also be used to put constraints on the total width of the Higgs boson [24,25]. In the light of the above motivations, we compute the full set of next-to-leading order (NLO) electroweak corrections to H → Z ( * ) Z ( * ) → 4 , with 4 = 4e, 4µ, 2e2µ. We match the NLO corrections to a QED Parton Shower (PS), in order to simulate multiple photon emission exclusively and provide final results at NLOPS electroweak accuracy. The calculation is available in an event generator, Hto4l 1 , which can be interfaced to any code describing Higgs boson production. The PS approach is based on the ideas first presented in Ref. [26] for the simulation of the Bhabha scattering process at GeV-scale e + e − colliders and later applied to the study of single W/Z production in hadronic collisions [27,28]. The matching procedure is a generalization of the method developed in Refs. [29,30] for the precision calculation of 2 → 2 processes in QED (as encoded in the program BabaYaga@NLO [31,32]) and also implemented in the event generator Horace for the calculation of single W/Z hadroproduction processes at NLOPS electroweak accuracy [33,34].
The NLO electroweak and QCD corrections to H → 4 fermions decay processes have been calculated in Refs. [35,36] and are available in the Monte Carlo (MC) program Prophecy4f [37,38], which is used in the context of Higgs studies at the LHC for the precision calculation of the branching ratios of the decays H → Z ( * ) Z ( * ) /W ( * ) W ( * ) → 4 fermions. In Prophecy4f higher-order photonic corrections are taken into account in terms of QED collinear Structure Functions. A preliminary study of the impact of the gauge-invariant NLO QED and PS corrections to the determination of the Higgs boson mass in the H → 4 decay was performed in Refs. [39,40].
The article is organized as follows. In Section 2 we describe the details of our calculation, with particular emphasis on the method used for the matching of the NLO electroweak corrections with the QED PS. In Section 3 we present our phenomenological results: in Section 3.1 we show a sample of comparisons between our predictions and those of Prophecy4f as a benchmark of the NLO computation, in Section 3.2 we provide results for various observables at NLOPS EW accuracy, while in Section 3.3 we present the results for Higgs production and decay in the channel gg → H → 4 obtained in terms of POWHEG [41,42] interfaced to Hto4l. In Section 4 we draw our conclusions.
2 Details of the calculation 2.1 Next-to-leading order (NLO) electroweak corrections The NLO electroweak corrections to the Higgs boson decay into four charged leptons consist of QED and purely weak contributions. Since the H → 4 decay is a neutral-current process, the two subsets are separately gauge invariant and can be computed separately as well.
The O(α) QED corrections are obtained by attaching a virtual or real photon to each charged lepton leg. They are expected a priori to provide the dominant contribution, as photons which are emitted collinear to a lepton give rise to large logarithmic corrections of the form α log m 2 /Q 2 , where m is the lepton mass and Q some typical energy scale.
The QED virtual corrections comprise vertex and pentagon diagrams (in the on-shell renormalization scheme), while real photon corrections are induced by the bremsstrahlung process H → 4 + γ. The two contributions are separately infrared (IR) divergent but their sum is IR-finite. We treat the IR singularity according to the standard QED procedure of assigning a small fictitious mass λ to the photon in the computation of the virtual and real contributions. More precisely, the Higgs decay width associated to the bremsstrahlung correction is separated in two pieces and calculated as follows (in a shorthand notation) where M H is the Higgs mass, is a soft-hard energy separator ( M H ), M 0 (H → 4 ) is the amplitude of the lowest-order (LO) process H → 4 and M(H → 4 + γ) is the matrix element of the radiative decay process H → 4 + γ, dΦ 5 being the 4 leptons plus 1 photon phase space element. In Eq. (2.1) eikonal factor stands for the analytical expression of the real radiation correction in the soft limit E γ → 0. The integral in the first line can be done analytically (see e.g. [43]) and the one in the second line is performed using standard MC integration with importance sampling.
The QED virtual counterpart is computed according to the following formula where M QED 1 is the one-loop amplitude associated to the O(α) vertex and pentagon diagrams.
We perform the IR cancellation by taking the sum of Eq. (2.1) and Eq. (2.2) in the numerical limit λ → 0. As a cross-check of the calculation, we tested that the inclusive NLO QED correction coincides with 2·3/4 (α/π), which is correctly twice the inclusive final-state O(α) electromagnetic correction to the Z → + − decay [44]. An important comment is in order here. The tree-level amplitude, as well as the amplitude for the real radiation process, contains poles in the phase space, corresponding to the points where the momenta of the + − pairs and of the ( + − γ) system cross the zero of the inverse propagators: These poles are avoided considering that the Z boson is an unstable particle, i.e. its propagator contains the finite Z-width. This, however, would spoil the IR cancellation between real and virtual corrections of Eq. (2.1) and Eq. (2.2), respectively, unless in Eq. (2.2) the QED virtual corrections are calculated with unstable Z bosons. The scheme which we adopt for the introduction of the width in the Z boson propagator, without introducing gauge invariance violations, is the complex mass scheme [45,46], which also allows us to include weak loop corrections consistently. 2 Concerning the basic features underlying the computation of the complete O(α) virtual corrections, we briefly describe the most important aspects in the following. Since we work in the 't Hooft-Feynman gauge, all the particles present in the spectrum of the SM, including the Fadeev-Popov and Higgs-Kibble ghosts, are involved in the calculation. The corresponding Feynman diagrams include, in addition to two-point functions, rank-two tensor three-, four-and five-point functions. The related ultraviolet divergencies are regularized by means of dimensional regularization. The reduction of the tensor n-point functions is carried out by means of the symbolic manipulation program FORM [47]. The necessary scalar form factors with complex masses are evaluated using Looptools v2.10 [48,49], which implements the evaluation of the reduction of tensor five-point integrals according to Refs. [50,51], as well as according to Passarino-Veltman reduction techniques [52]. The form factors are calculated with complex masses and real external squared momenta. This is sufficient for the implementation of the "simplified version of the complex renormalization", as described in Refs. [45,46]. The complete expressions for the counterterms in the on-shell scheme and for the basic self-energy diagrams are taken from Ref. [53]. Since the collinear singularities associated to the photon becoming collinear with one of the leptons are regulated by the finite lepton mass, the kinematics of the radiative process is calculated including exactly the contribution of lepton masses. In order to allow the cancellation of soft IR singularities, also the tree-level kinematics is calculated with complete lepton mass effects taken into account. In addition, this gives automatically the correct phase space integration boundaries for the diagrams of the virtual contribution where a virtual photon is connected to one external lepton pair. Although the kinematics is treated exactly, the non-IR O(α) virtual amplitudes are calculated in the approximation of neglecting finite fermion mass effects (with the exception of the quark Yukawa couplings, e.g. in the fermion-loop Higgs vertex corrections). These contributions are neglected in our calculation as they are irrelevant in view of a target theoretical accuracy of the order of 0.1% and their inclusion would make the numerical computation more time consuming.
In formulae, the Higgs width including one-loop weak corrections is obtained as where M weak 1 is the one-loop amplitude associated to the full set of O(α) weak diagrams. To check some relevant ingredients of our calculation of one-loop weak corrections, we compared our predictions for the Higgs decays H → ZZ, γγ at NLO electroweak accuracy with those of Ref. [54], finding perfect agreement.
In conclusion, our predictions for the Higgs boson decay into four leptons at NLO EW accuracy are given by the sum of Eq. (2.1), Eq. (2.2) and Eq. (2.3), supplemented with the necessary renormalization conditions.

Matching NLO electroweak corrections to QED Parton Shower
In the present section, we sketch our scheme for the matching of the NLO EW corrections with a QED PS. We closely follow the approach already presented and successfully applied to QED processes at low energies and Drell-Yan W/Z production at hadron colliders [29,30,33,34].
On general grounds, the partial decay width corrected for the emission of an arbitrary number of photons in a PS framework can be written as follows: where {p, k} stands for the set of the final state lepton and photon momenta p 1 , p 2 , p 3 , p 4 , k 1 , · · · , k n , |M P S n | 2 (of order α n ) is the PS approximation to the squared amplitude for the decay H → 4 + nγ, dΦ n is the exact phase space for the decay and Π({p}, ) is the Sudakov form factor accounting for unresolved emission, i.e. soft (up to a cut-off energy ) and virtual corrections in the PS approximation. It is understood that the integral over the phase space has a lower limit for the photon energies set to , to ensure the cancellation of the IR divergencies.
The quantities dΦ n and Π({p}, ) read explicitly In Eq. (2.6), L generates the soft/virtual collinear logarithms, including also interferences effects of radiation coming from different charged legs, and I , the integral of the Altarelli-Parisi vertex for the branching → +γ, generates the infrared logarithms. In the definition of L, the integral is performed over the angular variables of k, and η i equals 1 if i is an anti-fermion or −1 if it is a fermion. The integral over the phase space as in Eq. (2.5) is performed after choosing a convenient set of independent variables and using multi-channel MC importance sampling techniques to improve the integration convergence and follow the peaking structure of the partial decay width of Eq. (2.4) to help events generation. The fully exclusive information on final state particles momenta is kept. Details of the implementation are given in appendix A.
Before discussing the inclusion of NLO corrections into Eq. (2.4), it is interesting to point out that the squared amplitudes with photon emissions are enhanced in regions of the phase space where the photons are soft and/or collinear or where the Z propagators are resonating. In this perspective, a good approximation to the exact matrix element can be written in the form 3 : In the previous equation, c is a shorthand for the HZZ coupling, {σ, τ } label fermion and photon elicities, are the photon polarization vectors. P is a n-dimensional vector whose i th component is the index of the fermion to which the i th photon is attached and the sum over P denotes all possible ways to share n photons among the four fermions. Finally, Q P is the sum of the momenta of the photons, for a given P, attached to the electron current (R P to the muon current).
Equation (2.7) is derived from the amplitude for the emission of photons in the soft limit but keeping the dependence on the photon momenta in the Z propagators. The sum over the elicities of the squared amplitudes of Eq. (2.7) gives an approximation of the exact squared matrix elements, coherently including also interferences among diagrams. The final step to obtain |M P S n | 2 of Eq. (2.4) from Eq. (2.7) is to replace the photon energy spectrum with the Altarelli-Parisi distribution for a better treatment of hard collinear radiation.
Equation (2.4), with the building blocks described above, can then finally be improved to include exact NLO corrections according to our master formula: The (k i ) is its PS approximation. We want to remark that F SV and F H,i are by construction free of collinear and/or infrared logarithms and that the O(α) expansion of Eq. (2.8) exactly coincides with the NLO calculation, without any double counting. Furthermore, Eq. (2.8) is still fully differential in the final state momenta and can be conveniently implemented in a MC event generator.
Finally, we remark that the NLO virtual and real corrections used in F SV and F H,i are strictly defined only for 0 or 1 photon, while in Eq. (2.8) they are used also when there are additional photons: this requires a mapping of the n-photons phase space to 0 or 1 photon phase space. The mapping is implemented in close analogy to the one described in appendix A.2 of Ref. [29], and here we do not discuss it in further detail.

Numerical results
In the present Section, we show and discuss the numerical results provided by our calculation, as obtained with the new tool Hto4l. First, we show some tuned comparisons with the predictions of the reference code Prophecy4f at the level of NLO electroweak corrections. Then, we present our best predictions for various observables at NLOPS electroweak accuracy, as well as for Higgs production and decay in the presence of NLO QCD and electroweak corrections matched to PS.
The results presented here are obtained using Prophecy4f v2.0. 4 In both codes, we use the following set of input parameters  The M Z,W and Γ Z,W are the running-width PDG values which have to be converted to the fixed-width scheme adopted here through, for example, the relations of Eq. (7.2) of Ref. [35]. As we work in the G µ scheme, for the electromagnetic coupling constant we use the expression in the calculation of the LO width and NLO weak corrections, while we use α(0) for the coupling of the photon to the external charged particles. 5 The top-quark width is set to the LO prediction in the SM, and a fixed width is employed in all the resonant propagators in the framework of the complex mass scheme.

NLO electroweak corrections: comparisons to Prophecy4f
A sample of the Prophecy4f vs. Hto4l comparisons at NLO electroweak accuracy is shown in Tab. 2 and in Figs. 1-3, in order to check the technical accuracy of our predictions in its different aspects sketched in Sect. 2.1. Generally speaking, we observe very good agreement between our predictions and the independent results of Prophecy4f.
In Fig. 1   In Fig. 2 a comparison between Prophecy4f and Hto4l is shown for the e + e − invariant mass (in the Higgs rest frame), in the range [60,100] GeV (upper plot) and in the range [85, 95] GeV (lower plot). The results refer to the decay channel H → 2e2µ for M H = 125 GeV. Also in this case, the agreement between the two codes is remarkable, in spite of the large effect due to the radiative corrections 7 . Actually, at and above the peak of the electron-pair invariant mass distribution the corrections are of the order of 30%, while for M e + e − below M Z they can reach 50%. The lowering of the peak and the raising of a tail can 5 This value is used for all the numerical results shown in the following, with the exception of the comparisons with Prophecy4f, where we use αG µ everywhere, to be consistent with the default choice of Prophecy4f. 6 Analogous results are valid in the H → 4e channel, which coincides for the integrated partial width with the 4µ final state (apart from negligible mass effects). 7 For simplicity, here and in the following we provide results for bare electrons, i.e. in the absence of lepton-photon recombination effects. be mainly ascribed to the photon radiation off the leptons, as typical final-state radiation (FSR) effect observed around the peak of resonant processes [33,34,55,56].
A further comparison is given in Fig. 3 for the distribution of the angle between the decay planes of the virtual Z bosons in the H rest frame for the channels H → 2e2µ (upper plot) and H → 4µ (lower plot) for M H = 125 GeV, which is the observable of main interest for spin-parity assignment. For the φ angle we use the definition where k 12 = k 1 + k 2 and k 1 , k 2 , k 3 , k 4 are the three-momenta of the final-state leptons. Again the predictions of the two codes nicely agree. The contribution of the NLO corrections is particularly visible at the edges of the distribution, where it can reach the 5%  level for both the decay channels.

Predictions at NLOPS electroweak accuracy
Some illustrative results obtained according to a number of variants of the theoretical approach described in Sect. 2.2 are given in Figs. 4-5. In order to disentangle the impact of the different sources of correction, we consider the results obtained according to the following levels of accuracy: The comparison between approximations 1. and 2. is useful to quantify the higherorders contribution due to photon emission beyond O(α), while the difference between options 3. and 4. is a measure of pure weak loop corrections, the difference between approximations 2. and 3. is an estimate of non-logarithmic O(α) QED terms plus pure weak loop corrections. The comparison between approximations 3. and 5., as well as between 4. and 6., allows us to check that the NLOPS matching procedure correctly preserves the effect of QED exponentiation as given by the difference between options 1. and 2. Moreover, the results of 1. vs. those of 5. and of 3. vs. those of 5. provide an estimate of the accuracy of the predictions available in the literature for Higgs physics at the LHC, in particular of of the process-independent, widely used code PHOTOS [57], which describes multiple photon emission but does not include exact NLO electroweak corrections, and of Prophecy4f, that does not take into account the contribution of exclusive QED exponentiation.
In Fig. 4 we show the relative contribution of the different theoretical approximations discussed above for the e + e − (upper plot) and µ + µ − (lower plot) invariant mass in the Higgs rest frame, in the range [85, 95] GeV. The results refer to the process H → 2e2µ for M H = 125 GeV, according to a bare lepton definition. By inspection of Fig. 4 we can draw the following conclusions: the NLO corrections to the lepton invariant masses are quite large, since they amount to about 50% (30%) to the e + e − (µ + µ − ) invariant mass below the peak and about 30% (20%) at and above it. They are largely dominated by the enhanced leading logarithmic contributions of QED nature ∝ α log(M 2 Z /m 2 ), as can be inferred from the comparison between the results of the pure O(α) PS algorithm and those of the NLO QED/electroweak calculations. From this comparison, one can also conclude that the O(α) non-logarithmic QED terms contribute at the some per cent level, both for the e + e − and µ + µ − invariant mass, whereas the pure weak loops have a much smaller effect, not exceeding the 1% level.
The large impact of NLO QED corrections, which significantly modify the shape of the invariant mass distribution, translates in a relevant contribution due to higher-order photonic corrections. Multiple photon emission is of the order of 10% for the e + e − finalstate and at the level of some per cents for the µ + µ − case, as a consequence of the different  magnitude of the lepton-photon collinear logarithm. It can also be noticed that QED exponentiation reduces the impact of NLO corrections and that the NLOPS matching correctly preserves the size of multiple photon emission. Quite different conclusions derive from the analysis of Fig. 5, which shows the relative corrections of the different theoretical recipes on the φ angle distribution. For such an observable, the pure O(α) PS approximation significantly underestimates the contribution of NLO EW corrections for φ close to 0 • and 360 • , while it provides an overestimate around 180 • . Actually, it can be noticed that the φ angle distribution receives a non-  negligible contribution from fixed-order non-logarithmic terms and that, more importantly, is particularly sensitive to pure weak corrections, which set the correct overall size and shape of the radiative corrections. On the other hand, the effect of QED exponentiation is moderate, at the 1% level.
To summarize, the main conclusion is that both NLO electroweak and higher-order QED corrections, as well as their combination, are relevant for reliable simulations of the most important observables considered in precision studies of the Higgs sector at the LHC.

Interface to POWHEG: results for production and decay
In order to facilitate phenomenological studies of Higgs boson production and decay in the presence of both QCD and electroweak contributions, we have implemented an interface which allows to use our code in association with any event generator describing Higgs production. In Figs. 6-8 we show a sample of illustrative results obtained by interfacing Hto4l with POWHEG [42] for the simulation of Higgs boson production in gluon-gluon fusion. We use the POWHEG version with NLOPS accuracy in QCD [58] from the POWHEG BOX framework [59] and we consider Higgs production in proton-proton collisions at a c.m. energy of 8 TeV 8 . The events generated by POWHEG are interfaced to Hto4l according to the following procedure: • generate unweighted events for the process pp → H(+j) in the Les Houches format using POWHEG, where H is an on-shell Higgs boson and j stands for the extra parton of the NLO QCD calculation; • the Les Houches file is read event by event by Hto4l and the particles momenta are stored in the generic common block structure introduced in Ref. [60]; • each event is decayed into the selected channel in the H rest frame, using Hto4l.
After boosting the decay products back to the laboratory frame, the events including production and decay are written in a file in the Les Houches format.
The Les Houches file can be finally passed to a shower event generator for QCD showering and hadronization. In our examples we use PYTHIA v6.4 [61] as QCD PS. According to the above procedure, the pp → H → 4 process is treated in narrow width approximation, as it is the case for a 125 GeV Higgs boson, and factorized in on-shell Higgs production and decay.
In our analysis we consider, for definiteness, the decay channel H → 2e2µ and the following observables: the transverse momentum p H T and rapidity y H of the Higgs boson (Fig. 6), the invariant mass of the subleading lepton pairs and the magnitude of the cosine of the decay angle of the leading lepton pair in the four-lepton rest frame with respect to the beam axis | cos θ * | (Fig. 7). The leading pair is defined as the lepton pair with invariant mass closest to the Z boson mass and its angle is obtained by summing the three-momenta of the two leptons. For the POWHEG calculation of Higgs production in gluon fusion, we use the PDF set MSTW2008nlo68cl [62] with factorization/renormalization scale µ R = µ F = M H . The values of the other input parameters are the same as the ones given in Tab. 1. The results shown in the following refer to a sample of 1, 2 · 10 8 unweighted events and to the same selection cuts adopted in Ref. [23] and correspond to bare-level leptons.
In Fig. 6 and Fig. 7 we show the comparison between the predictions obtained using POWHEG interfaced to our code at LO and NLOPS electroweak accuracy. It can be noticed that the contribution due to NLOPS electroweak corrections is almost flat and of about −15% for p H T , y H and | cos θ * |, while the invariant mass of the subleading lepton pairs receives a varying correction of size between −20% and −10%.
In Fig. 8 we show the results for two observables which are fully exclusive over QED radiation and which can be easily treated in our approach. The results correspond to the process pp → H → 2e2µ + n γ, with E min γ = 1 GeV, for which we show the transverse momentum of the hardest photon and the angular separation between the hardest photon and the closest lepton, that exhibit the expected features of photon emission in radiative events.

Conclusions
In this work we have presented a precision calculation of the SM Higgs boson decay into four charged leptons, in view of improved measurements of the properties of the Higgs particle at the LHC Run II. Our approach is based on the computation of the full one-loop electroweak corrections supplemented with the contribution of multiple photon emission taken into account according to a fully exclusive QED PS algorithm. Our results, which have a NLOPS electroweak accuracy, are available in the form of a new event generator, Hto4l, that can be easily interfaced to any QCD program simulating Higgs production. We have cross-checked our NLO electroweak corrections against the predictions of the reference code Prophecy4f and found perfect agreement. We have also shown that both NLO electroweak and higher-order QED corrections, as well as their interplay, are necessary for actually precise simulations of the variety of observables involved in Higgs physics at the LHC. This provides the main novel theoretical feature of our work, which goes beyond the presently available results limited to the fixed-order approximation or to a leading logarithmic QED modeling. The second relevant aspect is given by the possibility of interfacing Hto4l to any generator describing Higgs boson production, thus allowing simulations of Higgs production and decay in the presence of higher-order QCD and electroweak correc-tions. In this respect, we have shown some illustrative results obtained in terms of the combined usage of POWHEG and Hto4l.
Our results can find application in precision measurements of the Higgs boson mass, spin-parity determination and tests of the SM at the level of differential cross sections in the future run of the LHC. They can be generalized to other processes yielding four leptons in hadronic collisions, like e.g. pp → H → W ( * ) W ( * ) → 2 2ν or pp → Z ( * ) Z ( * ) → 4 .

A Phase space parameterisation and integration
The 4 + n bodies phase space as in Eq. (2.5) is integrated according to standard multichannel MC techniques, combined with importance sampling to reduce the variance of the integral and help event generation 9 . The first step is to generate a photon multiplicity n and associate n 1 (n 2 ) photons to the electron (muon) current (n 1 + n 2 = n), defining the channel of the multi-channel integration. The phase space is then conveniently split into two decaying objects to follow the Z propagators, namely dΦ(P H ; p 1 , · · · , p 4 , k 1 , · · · , k n ) = (2π) 6 dQ 2 Z 1 dQ 2 Z 2 dΦ(P H ; P Z 1 , P Z 2 ) × dΦ(P Z 1 ; p 1 , p 2 , k 1 , · · · , k n 1 ) dΦ(P Z 2 ; p 3 , p 4 , k n 1 +1 , · · · , k n 1 +n 2 ) (A.1) where P Z i (P 2 Z i = Q 2 Z i ) are the momenta of the virtual Z bosons. We refrain from writing explicitly the simple 1 → 2 decay phase spaces of Eq. (A.1) and we focus instead on the case where at least one photon is present. As discussed in appendix A.3 of Ref. [29], an efficient sampling of photons collinear to final state leptons is a non trivial task, because the directions of the leptons are known only after all the momenta are generated. In Ref. [29] we adopted a solution based on a properly chosen multi-channel strategy. Here we adopt a different and elegant solution, which consists in writing the phase space in the frame where the leptons are back-to-back, i.e. p a = − p b (see for example [63][64][65]).
Omitting overall numerical factors for brevity, the building block we are interested in is dΦ(P ; p a , p b , k 1 , · · · , k r ) = δ (4) Here we consider only the decay H → 2e2µ, the generalization to 4 identical leptons being straightforward.
where we defined Q = p a + p b , K = r i=1 k i and δΦ contains the infinitesimal phase space element divided by the final state particle energies. It is usually understood that all the variables are expressed in the frame where P is at rest, but we want to express them where Q is at rest. In order to do that, the previous equation can be further manipulated by inserting the following identities which help to make explicit the Lorentz invariance of the phase space element.
With the help of Eq. (A.2) and appropriately rearranging the terms, we can write dΦ(P ; p a , p b , k 1 , · · · , k r ) = δΦ d 3 Qδ (4) In the cascade of identities (A.3) we used the result d 3 Q δ (3) ( P ) = (s /s) 3 2 (see [64]) and we made use of Lorentz invariance. In the last identity it is understood that all the variables are expressed in the frame where Q = p a + p b is at rest and s = Q 2 , s = P 2 , β a is the speed of particle a and dΩ a = d cos θ a dφ a . The big advantage of the last equation is that the lepton momenta p a and p b lie on the same direction defined by cos θ a and φ a , hence all photons can be generated along this direction to sample the collinear singularities. Once all particle momenta are generated, they can be boosted back to the rest frame of the decaying Higgs boson.
One last remark concerns the integration limits of the phase space. As mentioned in Sect. 2.2, photon energies should be generated larger than the infrared cut-off in the Higgs frame, which is a non Lorentz invariant cut. Since the minimum photon energy can not be determined a priori in the frame where Q is at rest (because Q itself depends on the photons momenta), we decide to generate photon energies starting from 0 to cover the whole phase space and then, once boosted back, cut the event if a photon enegy falls below . Finally, in order to flatten the infrared divergence, we choose to sample the photon energies according to the function where is a guessed (and tuned for efficiency) minimum energy.