Higgsstrahlung at NNLL$'+$NNLO Matched to Parton Showers in GENEVA

We present results for the Higgsstrahlung process within the GENEVA Monte Carlo framework. We combine the fully differential NNLO calculation with the higher-order resummation in the 0-jettiness resolution variable (beam-thrust). The resulting parton-level events are further showered and hadronised by Pythia8. The beam-thrust resummation is carried out to NNLL$'$ accuracy, which consistently incorporates all singular virtual and real NNLO corrections. It thus provides a natural perturbative connection between the NNLO calculation and the parton shower regime, including a systematic assessment of perturbative uncertainties. In this way, observables which are inclusive over the additional radiation are correct to NNLO, while the description of 0-jet-like resummation variables is improved beyond the parton shower approximation. We provide predictions for the 13 TeV LHC.


INTRODUCTION
The study of the properties of the Higgs boson forms an important part of the experimental programme at the LHC. Though first observed in final states with two photons, the dominant decay mode of the Higgs boson is in fact to a pair of bottom quarks which has only recently been measured in associated production with a vector boson [1,2]. This class of Higgsstrahlung processes is, therefore, an interesting subject of study. In particular, it allows one to probe couplings of the Higgs boson to electroweak vector bosons as well as to heavy quarks. It is also a particularly clean channel, since the leptonic signature arising from the bosonic decay can be efficiently distinguished above the QCD background.
At tree level the calculation involves a relatively simple extension of the Drell-Yan process in which a Higgs boson is emitted from the vector boson. In the case of ZH production, the gluon-initiated channel in which the Higgs boson couples to a top-quark loop becomes accessible beginning at O(α 2 s ). This provides a finite contribution to the next-to-next-to-leading order (NNLO) cross section and is particularly sensitive to the presence of New Physics. Other top-quark loop mediated contributions, while in principle present, have been shown to provide only a percent level contribution to the total cross section [3] and are often neglected in differential calculations, being separately finite and gauge invariant.
There has also been considerable progress in the development of Monte Carlo event generators. NLO calculations matched to parton showers for Higgsstrahlung have been available in the Powheg framework [19] for quite some time. More recently merged NLO calculations of V H and V H+jet have appeared [19], utilising the MiNLO method [20,21]. Of the available methods which can attain NNLO accuracy matched to parton shower [21][22][23][24], to date only the MiNLO approach has been applied to the Higgsstrahlung case [25,26]. Merged NLO samples of V H and V H+jet including NLO EW corrections matched to parton showers also produced using the MiNLO method have been presented in Ref. [27].
In this paper, we present results for the Higgsstrahlung processes pp → ZH → + − H and pp → W ± H → ± ν H at NNLO matched to the Pythia 8 parton shower according to the Geneva method [23,[28][29][30], assuming an on-shell Higgs boson. The fully differential Higgsstrahlung calculation at NNLO is improved with next-to-next-to-leading logarithmic (NNLL ) resummation of 0-jettiness (beam thrust) and subsequently showered while maintaining NNLO accuracy at FO for the underlying process.
The content of the paper is organised as follows. In sec. 2 we present a review of the Geneva framework in some detail, in order to provide background for the uninitiated reader. We discuss the implementation of the Higgsstrahlung processes and highlight details specific to this case. Next, in sec. 3, we describe the validation of our results against existing NNLO predictions. In sec. 4 we then present our predictions for various distributions at partonic, showered and hadronised levels. Finally, we make some general comments on the outcome of the work and present potential future directions in sec. 5.

THEORETICAL FRAMEWORK
The theoretical framework of our calculation is based on the Geneva method, which has been developed in Refs. [23,[28][29][30]. While we report in the following all the necessary ingredients for the specific Higgsstrahlung process under investigation, we refer the interested readers to the original papers where the derivations of the results used here are presented in greater detail.

A. General Setup
The Geneva framework is based on the definition of physical and infrared (IR)-finite events, generated at a given perturbative accuracy and obtained from both FO and resummed calculations. This is achieved by translating IR-divergent final states with M partons into IRfinite final states with N partonic jets (where M ≥ N ) for which the divergences cancel on an event-by-event basis. The translation is performed using an N -jet resolution variable T N which partitions the phase space into regions with different numbers of resolved emissions in the final state. For example, the Geneva 0-jet cross section dσ mc 0 receives contributions not only from 0-parton events, but also from 1-parton events where the additional emission is unresolved, i.e. below the T cut 0 value of the 0-jet resolution, and from 2-parton events where both additional emissions are unresolved. The partitioning of the phase space is achieved by defining cuts on the T 0 and T 1 resolution parameters and separating regimes as follows: 1 In this way the cross section over the entire phase space is divided into exclusive 0-jet and 1-jet cross sections and an inclusive 2-jet cross section. The partonic jet bins thus defined are rather different from those an experimentalist might define using a usual jet algorithm: their definition depends on an IR-safe phase space map Φ N (Φ M ) which projects an M -body onto an N -body phase space and ensures the individual IR-finiteness of the resulting MC cross sections dσ mc i . Using the events in eq. (1), the total cross section for an observable X is then given by where M X (Φ N ) is the measurement function that computes the observable X for the N -parton final state Φ N . 1 We exploit the notation in eq. (1) to highlight the dependence of the dσ mc i cross section on the resolution parameters. When an argument contains a single term, e.g. T cut N , it means that the corresponding quantity has been integrated over up to the value of the argument. An argument T N > T cut N implies instead that the corresponding cross section is still differential in the relevant resolution variable for values larger than the cutoff. Despite appearances, this object is not identical to the result one would obtain in a standard FO calculation since the unresolved emissions are assigned to the projected phase space points Φ N (Φ M ) rather than to the exact points Φ M . However, provided that the definitions of the dσ mc i cross sections correctly capture all the singular contributions at the given order, the nonsingular difference due to the projection vanishes in the limit T cut N → 0. We are therefore motivated to choose as small a value for T cut N as possible. We should then expect to encounter large logarithms of T N and T cut N as we take this limit, and, to prevent the convergence of our perturbation theory being spoiled, these terms must be resummed.
The attentive reader may have noticed that the requirement of projectability from Φ M → Φ N allows for nonsingular (non-projectable) events to be assigned to the higher multiplicity jet bins despite their value for the T N resolution variable being below the T cut N cutoff. A simple example can be seen in Φ 1 configurations qg → ZHq: when the direction of the outgoing quark is collinear to that of the incoming quark (and therefore anti-collinear to the incoming gluon) the resulting Φ 0 projection would require a gg → ZH tree-level configuration which does not exist. These events are therefore classified as Φ 1 and, lacking any collinear or soft enhancement due to their nonsingular nature, they are assigned the corresponding tree-level cross section.
The separation between the 0-jet and 1-jet regimes is determined by a 0-jet resolution variable T 0 . In Geneva, resummation of this variable is performed at NNLL accuracy, therefore including all contributions singular in T 0 to O(α 2 s ). 2 The exclusive 0-jet and inclusive 1-jet cross sections can then be written as where the nonsingular terms contain at worst integrable singularities. In the case of the inclusive 1-jet cross section, it is necessary to extend the dependence of dσ NNLL from Φ 0 to Φ 1 by including the differential dependence on the radiation phase space, parameterised in terms of T 0 and two other variables. This has been done by considering the resummed differential spectrum in T 0 , dσ NNLL /dT 0 dΦ 0 , and introducing a normalised splitting function P(Φ 1 ) to account for the dependence on the two remaining variables. These could be for example the fractional energy z of one daughter in the splitting and an azimuthal angle ϕ. In order not to spoil the normalisation for each point in the T 0 spectrum, the splitting function must satisfy Since we wish to obtain overall NNLO accuracy, we must have dσ mc 0 and dσ mc ≥1 at NNLO 0 and NLO 1 , respectively, which determines the nonsingular matching contributions to be The terms in square brackets are the FO expansions to O(α 2 s ) of the resummed cumulant and spectrum. Inserting the expressions for the FO cross sections, we obtain dσ mc where B M contains the M -parton tree-level contributions, V M the M -parton one-loop contributions, and W 0 the two-loop contribution. We have also introduced the shorthand notation Because the resummed contribution is differential in T 0 , particular care is needed when integrating the FO 2parton contribution B 2 of the inclusive 1-jet cross section in eq. (9). The associated radiation phase space must be parameterised specifically by T 0 and two other arbitrary variables, e.g. dΦ 1 = dΦ 0 dT 0 dz 1 dz 2 where the z i might for example be {z, ϕ}. The projection dΦ 2 /dΦ T 1 , which implicitly defines Φ T 1 , must therefore use a map which preserves the value of T 0 , such that the pointwise singular T 0 dependence is alike among all terms in eq. (9) and cancellation of these singular terms is guaranteed. The projection used is defined as where Θ T (Φ 2 ) defines the region of Φ 2 that can be projected onto the physical Φ 1 phase space via the IR-safe map Φ T 1 (Φ 2 ). Only this projectable region of Φ 2 is included in dσ nons ≥1 /dΦ 1 , while the remainder will be included in the nonsingular Φ 2 events below.

1/2-jet separation
The separation of the inclusive 1-jet cross section into an exclusive 1-jet cross section and an inclusive 2-jet cross section proceeds in analogy to the 0/1-jet case, with the relevant resolution variable now T 1 and the requirement on the resummation accuracy relaxed to NLL. We write dσ mc where now, in contrast to the 0/1-jet case, it is sufficient to consider contributions only up to NLL to ensure that the matching terms, which determine the FO accuracy, are free from singular logarithmic enhancements in T 1 . For NNLO accuracy, dσ mc 1 and dσ mc ≥2 must be correct to NLO 1 and LO 2 respectively, and so the nonsingular matching contributions are (c.f. eqs. (6) and (7)) At NLL, the resummed contributions take the form We now discuss separately each contribution appearing in the equations eqs. (17) and (18) above. In order to understand the rôle of dσ C ≥1 , it is perhaps clearer if we temporarily neglect the T 0 resummation which was constructed at NNLL accuracy in the previous subsection. In this approximation, dσ C ≥1 is a proxy for dσ FO

≥1
in the limit that T 1 → 0. We choose it to have a similar form to the FO contribution to eq. (9) but with the full double-real matrix element B 2 replaced with its singular approximant C 2 (Φ 2 ) and with a different projection, . The expanded form of dσ C ≥1 at NLO 1 is therefore given by The C 2 term acts as a standard NLO subtraction that reproduces the pointwise singular behaviour of B 2 -in practice, we have implemented the FKS subtractions [31].
U 1 (Φ 1 , T cut 1 ) denotes the Sudakov factor that resums the dependence on T cut 1 to NLL accuracy and U 1 (Φ 1 , T 1 ) denotes its derivative with respect to T cut 1 . At this order, it is given by [32,33] with Γ the Euler gamma function and The functions appearing in the above are defined as with r = α s (µ 2 )/α s (µ 1 ) and the dependence on T cut 1 appears via the dependence on the scales Here, T max 1 is the value at which the T 1 resummation is turned off, which is chosen near the maximum kinematically allowed value of T 1 for a given phase space point Φ 1 . The cusp and non-cusp anomalous dimensions entering the above expressions are well known, T F n f , while the kinematical terms are where p a , p b , and p 1 are the massless four-momenta of the Φ 1 phase space point, and The assignment of p a , p b , and p 1 to p q , pq, and p g is according to the flavour structure of Φ 1 . For example, for a qq → HZg flavour structure we have p q = p a , pq = p b and p g = p 1 . Inserting now eqs. (17) and (19) into eq. (15), we find that the matching term for the exclusive 1-jet cross section is given by while analogously in the 2-jet case we obtain In these expressions, U are the O(α s ) terms in the expansions of the objects U 1 and U 1 , which cancel the logarithmic terms in T 1 in the B 2 pieces. It is worth noticing that the contributions in eqs. (27) and (28) are actually nonsingular at O(α s ), despite the accuracy of the T 1 resummation being only NLL and not NLL . This is due to the fact that we have included the full O(α s ) virtual, soft and collinear contributions in eq. (19).
Thus far, we have discussed the construction of an additive NLO 1 +NLL T1 matching but have neglected to include the T 0 resummation which we constructed at NNLL in the previous subsection. In order to include this correctly we must ensure that the integral of the dσ mc 1 and dσ mc ≥2 cross section reproduces the T 0 -resummed result for the inclusive 1-jet MC cross section dσ mc ≥1 in eq. (9). That is, where we have used the identity (noting that U 1 (Φ 1 , T max 1 ) ≡ 1 and c.f. eq. (5)) Inserting the expression for dσ mc ≥1 in eq. (9) into eq. (29), we obtain the result for dσ C ≥1 beyond NLO 1 and thus the full expressions for the exclusive 1-jet and inclusive 2-jet cross sections, These contain the differential T 0 resummation via dσ NNLL and completely define the fully differential jet cross sections.
B. Implementation of the V H process in GENEVA

Choice of the jet resolution variables
We use N -jettiness [32] as our N -jet resolution variable, defined as where the sum over k runs over all coloured final-state particles and whereq i = n i = (1, n i ) are light-like reference vectors along the jet and beam directions. While the reference vectors which lie along the beam directions are the same for any N so that we can choose n a =ẑ and n b = −ẑ, for values of N ≥ 1 the definition of the reference vector along the jet direction depends on a clustering metric. We refer the interested reader to [23] for details.
N -jettiness quantifies the degree to which the final state is N -jet-like for a given N , and has the useful property that T N = 0 in the limit that a configuration is composed of exactly N partons. It can be used to cluster the final state into N -jet and beam regions in an IR-safe manner without resorting to any additional clustering algorithms. Crucially for our purposes, both its singular structure and resummation are known to the requisite accuracy.
For production of any colour singlet at NNLL , the two resolution variables 0-jettiness T 0 and 1-jettiness T 1 are needed to partition the phase space. Since this construction is identical to what has been previously done for the Drell-Yan process [23] we refer the reader to the discussion therein.

The T0 spectrum at NNLL from SCET
The all-orders parton-level factorisation theorem for 0jettiness is given by [34,35] where dσ B ij /dΦ 0 is the Born cross section for the The hard function H ij (Q) contains the corresponding Born and virtual squared matrix elements, and the sum runs over all possible qq pairs ij = {uū,ūu, dd,dd, . . .}.
The B i (t, x) are inclusive (anti)quark beam functions [34], with partonic virtualities t a,b and momentum fractions x a,b given in terms of the total rapidity Y V H and invariant mass Q = M V H of the V H final state as well as the hadronic centre-of-mass energy E cm by The beam functions are computed perturbatively in terms of standard PDFs f j ; schematically B i = j I ij ⊗ f j where the I ij are perturbative coefficients. Finally, S(k) is the quark hemisphere soft function for beam thrust [36].
The preceding eq. (35) is derived using Soft Collinear Effective Theory (SCET) [37][38][39][40][41][42][43]. We note that each of H, B, S depends only on a single characteristic scale. The consequence is that the perturbative expansions of these constituent functions do not feature any large logarithms when a suitable choice of scale is made, viz.
In eq. (35), however, all ingredients must be evaluated at an arbitrary common scale µ, whose dependence exactly cancels between the different functions at the appropriate order. We achieve this by using the renormalisation group evolution in the effective theory to evolve each function from its own scale to µ. We thus obtain the resummed T 0 spectrum used in eq. (7): where the large logarithmic terms arising from the ratios of scales have been resummed by the renormalisation group evolution (RGE) factors U X (µ X , µ). At NNLL accuracy, the boundary conditions for the evolution of the functions are required at 2-loop order, while the evolution kernel itself must be inserted at 3(2)-loop order in the cusp (non-cusp) anomalous dimensions. It suffices to say that all required expressions are known [44][45][46][47][48][49][50][51][52][53]; they are in fact mostly identical to those used for Drell-Yan production.
The only exception is the hard function, for which many of the relevant Feynman diagrams are still closely related to those appearing in the Drell-Yan case. In particular, for the class of diagrams illustrated in fig. 1a the loop corrections are identical and differences are due solely to the emission of the Higgs boson from the finalstate vector boson. Indeed, writing the hard function as the product of a hadronic and an electroweak tensor, it is apparent that the hadronic tensor H µ ij (where the indices i, j run over the allowed partonic flavours) is identical in the Higgsstrahlung and Drell-Yan cases for the class of diagrams in fig. 1a and only the electroweak tensor is modified. Since in Ref. [23] the hard function is implemented as a factor multiplying the squared Born amplitude, we may use that same factor here and simply replace the tree level piece with the appropriate Higgsstrahlung contribution. In the case of ZH production at NNLO, however, a second class of diagrams appears, depicted in fig. 1b. These gluon-initiated contributions are finite and enter at O(α 2 s ); due to the dominance of the gluon PDF at the LHC it is important that they are included, as they have an O(10%) effect on the total cross section [26]. As their contribution is purely nonsingular in nature, they do not affect the implementation of the resummation and can be included separately.
In the present calculation we neglect 2-loop contributions involving top quarks (which are separately finite and gauge-invariant) as their effect has been shown to contribute to the total cross section only at O(1%) [3] and their exact form remains unknown. However, we include real-virtual corrections in which the Higgs boson couples to a top-quark loop with exact top-quark mass dependence.

Scale choices and profile scales
The purpose of the resummation is to correctly account for the effects of the logarithms of lar and nonsingular terms, preventing one from obtaining the correct FO result. It is therefore important to switch off the resummation in the calculation before this occurs. In order to determine the relevant value of T 0 , it is instructive to consider the relative sizes of the singular and nonsingular contributions as functions of T 0 . We show the absolute values of these contributions in fig. 2. We see that these become similar in magnitude at T 0 ≈ 150 GeV and therefore turn off the resummation around this point. While this consideration holds for the total spectrum, integrated over all the possible Born-like kinematics, one may wonder exactly what the kinematical dependence on the coefficients of the logarithms is and whether this dependence might affect where the resummation needs to be turned off. We consider this possibility in fig. 3, where we plot the singular and nonsingular spectra as functions of τ 0 = T 0 /M V H . The figure shows that the location of the crossing point between the spectra exhibits modest dependence on the exact kinematic regime. We may therefore switch off the resummation in a similar manner for all classes of event.
Since the resummation is achieved via RGE running, it is sufficient to set all resummation scales to a common nonsingular scale, viz. µ NS = µ S = µ B = µ H , to stop the evolution. In order to ensure a natural transition between the resummation and FO regimes, we make use of profile scales µ B (T 0 ) and µ S (T 0 ) which are constructed to interpolate smoothly from the characteristic scales to µ NS [54,55]. We have where the common profile function f run (x) is as in Ref. [56], This form has strict canonical scaling (cf. eq. (37)) below x 1 and switches off the resummation above x 3 . Considering fig. 2, we are led to the choice of parameters In the resummation region, the nonsingular scale µ NS must be chosen of the same order as the hard scale Q for the inclusive Higgsstrahlung. In the FO region, it can instead be left free to match any arbitrary fixed or dynamic scale µ FO . The transition is achieved by imposing that Q = M V H for values of τ 0 up to x 3 and smoothly interpolating the µ NS value between Q and µ FO above that threshold.
We may estimate the uncertainties associated with the resummed and FO calculations by varying the profile scales. In the FO case, we adopt the usual prescription of varying µ NS up and down by a factor of 2 and taking the maximal absolute deviation from the central value as a measure of the uncertainty. This preserves everywhere the ratios between the various scales µ H , µ B and µ S and so the arguments of the logarithms which are resummed by the RGE factors are unaffected. In the resummed case, we vary independently the profile scales for µ B and µ S about their central profiles while keeping µ H = µ NS fixed. In this way the arguments of the resummed logarithms are varied in order to estimate the size of higher-order corrections in the resummed series while maintaining the scale hierarchy µ NS ∼ µ H µ B ∼ √ µ H µ S µ S . More details on the specifics of this prescription may be found in Ref. [57]. In addition, we include two more profiles where we vary all x i transition points by ±0.05 simultaneously while keeping all the scales at their central values. We thus obtain six profile variations in total and take the maximal absolute deviation in the result from the central value as the resummation uncertainty. The total perturbative uncertainty is then obtained as the quadrature sum of the resummation and FO uncertainties.
Hitherto we have considered only the effect of the scale choice on the T 0 spectrum, dσ NNLL /dT 0 , and not the effect on its integral over T 0 , the cumulant dσ NNLL (T cut 0 ). Indeed, while summing over the events distributed according to the T 0 spectrum, for example for the calculation of the inclusive total cross section, one performs exactly this integration. However, since the profile scales have a functional dependence on T 0 , the integral of the spectrum is not exactly equal to the cumulant evaluated one would obtain where T max 0 is the upper kinematical limit, so that the difference is due to terms of higher order. However, these terms could be numerically relevant, especially if one aims at reproducing the exact FO inclusive cross section.
In order to remedy this, we may add a term to the spectrum such that the inclusive FO cross section is correctly recovered. The term added must satisfy the following properties: • The integral of the modified spectrum must recover the FO cross section; • The term must contribute only in the region of T 0 where the missing N 3 LL terms are sizeable and must vanish elsewhere (especially in the FO region at large T 0 ); • The term must be of the same order as the missing terms so that the NNLL accuracy of the spectrum is not spoiled when the term is added.
We therefore add the term where κ(T 0 ) and µ h (T 0 ) are smooth functions. First, we note that the term is, by construction, of higher order and thus our third criterion is satisfied. Second, we note that in the FO region where µ h (T 0 ) = Q the difference between the pieces in brackets is zero, as the scales are constant here -the term therefore vanishes. We may also choose κ(T 0 ) to tend to zero in this region to reduce further the size of this additional contribution before exact cancellation is reached, and choose the profile scale µ h (T 0 ) to reach Q at a lower value of T 0 than the rest of the calculation. All this ensures that the accuracy of the tail of the spectrum is not spoiled by the addition of eq. (45) and instead restricts it to act in the region where T 0 ∼ T cut 0 Q. Our second criterion is therefore satisfied. Third, we may tune κ(T 0 → 0) such that, upon integration of the sum of the spectrum and eq. (45), the correct inclusive cross section is recovered. We thus satisfy our first criterion. In fact, we can perform this tuning for each FO scale variation separately. When taking either the central value of µ(T 0 ) or any of the resummation variations, we simply take µ h (T 0 ) as before. When we take the FO up and down variations for µ(T 0 ), however, we now take µ up h (T 0 ) = 2µ h (T 0 ) or µ down h (T 0 ) = 1/2µ h (T 0 ) and readjust the value of κ(T 0 → 0) such that the value of the inclusive cross section is correctly recovered. In this way we obtain the correct FO scale variations for the inclusive cross section (or indeed any other inclusive quantity).

Power-suppressed corrections to the nonsingular cumulant
The 0-jet cross section in eq. (8) is NNLO accurate and fully differential in the Φ 0 phase space. In order to regulate the IR divergences that appear in 4 dimensions in the intermediate stages of the calculation, however, one needs a subtraction method up to NNLO which is also fully differential in Φ 0 . While fully general local subtractions are available at NLO [31,58,59], local NNLO general subtraction methods are still in their infancy [10,[60][61][62][63][64]. Even when provided with a NNLO subtraction, the Geneva predictions at O(α 2 s ) stemming from eq. (8) would be exact for the total cross section but might be correct only up to power corrections in T cut 0 for observables dependent on the Φ 0 kinematics. This is a consequence of the projective map used in the definition of the Φ 0 events. 3 For this reason, one can avoid the necessity of a NNLO subtraction which is fully differential in Φ 0 and replace the formula for the 0-jet cross section with dσ mc where only a local NLO subtraction and the expansion of the resummation at O(α s ) are now required. The formula above assumes that there is an exact cancellation of both the singular and nonsingular contributions at O(α 2 s ) between the FO and the resummed-expanded terms. In reality these terms differ by a nonsingular contribution, which can be written as where the functions f i (T cut 0 , Φ 0 ) are at worst logarithmically divergent in the small T cut 0 limit. While we include the NLO term proportional to f 1 (T cut 0 , Φ 0 ) in eq. (46), we neglect the f 2 (T cut 0 , Φ 0 ) piece. The size of this neglected term as a function of the cut is shown in fig. 4a for ZH production. We note that at the default value at which we run, T cut 0 = 1 GeV, the observed effect is rather small in magnitude at just under 0.2 fb (∼ 0.7% of the total cross section). In order to correct for this discrepancy and obtain the correct NNLO inclusive cross section, we may simply rescale the weights of the Φ 0 events in such a way that the total cross section thus obtained matches the result of an independent NNLO calculation. In practice, we use cross sections obtained from the Matrix program [65] to obtain the reweighting factors for the central scale and its FO variations. By following such a procedure, we are able to include the effects of the f 2 term in eq. (47) on the total cross section that would have been present had we implemented eq. (8) literally. Since neither eq. (8) nor our approach in eq. (46) achieve the exact O(α 2 s ) Φ 0 dependence of all observables, our approximation does not inherently limit the accuracy of our predictions.
In nearly all spectra, the lack of the correct O(α 2 s ) Φ 0 dependence does not produce striking differences when compared with Matrix. We therefore conclude that our approximation holds rather well. As an example, we show the transverse-momentum distribution of the hardest lepton produced for different values of T cut 0 in fig. 4b. We observe similar behaviour in many other distributions. In one exceptional case, however, we find a mild effect on a distribution, namely the Higgs boson transverse momentum. In fig. 4c, one can see a difference of O(1%) between the Geneva and Matrix results in the first few bins. This difference is halved when the T cut 0 is reduced to 0.1 GeV and is restricted in range compared to the results using the higher cut. Nonetheless, throughout this work we have continued to use a value T cut 0 = 1 GeV for reasons of improved numerical stability (the reader may note the increase in the sizes of the statistical errors associated with the T cut 0 = 0.1 GeV calculation, for example). We find that any discrepancies caused by this choice are consistently small and would most likely fall within PDF uncertainties.

NLO1 calculation and phase space mapping
We remark that the approximation described above is applied only to the 0-jet bin. The Φ 1 and Φ 2 events are produced instead following the full forms in eqs. (32) and (33), which employ a proper 0-jettiness subtraction at NNLO combined with the local NLO 1 subtraction for the V Hj process. In practice we use the FKS subtraction, but with a specific choice for the mapping employed in splittings and projections. Since the FO terms appearing in eqs. (32) and (33) are differential in T 0 , the mapping used in the NLO 1 calculation must preserve the value of T 0 during such splittings or projections, see e.g. eq. (12). Such a mapping was devised for the Drell-Yan production in Ref. [23] and is discussed briefly therein. For this work, we performed only minor modifications in order to accommodate the additional final-state Higgs boson.

Interface to the parton shower
Since the shower interface is mostly identical for Drell-Yan production and Higgsstrahlung, we provide only a brief recap of its main features here, referring the interested reader to section 3 of Ref. [23] for a more detailed discussion.
The partonic jet cross sections dσ mc 0 , dσ mc 1 and dσ mc ≥2 each include contributions from higher multiplicity phase space points, but only in those cases where T N (Φ M ) < T cut N . In order to make the calculation fully differential in the higher multiplicities, a parton shower is interfaced which adds radiation to each jet cross section in a manner which does not alter the value or accuracy of the integrated cross sections but still produces a fully exclusive final state, i.e. in a unitary and recursive fashion. The effect of this is to restore the emissions in dσ mc 0 and dσ mc 1 which were integrated over when the jet cross sections were constructed, as well as to add extra final-state partons to the inclusive dσ mc ≥2 . In order to simplify the treatment, we consider a shower strongly ordered in T N , such that . . . In practice, no current parton shower program uses T N as an evolution variable in the way that we have here, choosing instead e.g. the transverse momentum of an emission. However, if one were to take the output of a shower ordered in say transverse momentum, one could recluster the partons using the Njettiness metric in order to obtain a splitting history that was ordered in T N and equivalent at LL order.
The requirement of the preservation of the accuracy of the jet cross section after applying the shower on a phase space point Φ N sets constraints on the point Φ N +1 reached after each emission added by the shower. These constraints are different for the different partonic multiplicities of the events before the shower.
For the cases in which the showered events originate from Φ 0 events, the main constraint is that the integral of the cross section below the T cut 0 (which is NNLL +NNLO accurate) must not be modified. The emissions generated by the shower must in this case satisfy T 0 (Φ N ) < T cut 0 , so that they recover the events which were integrated over in the construction of the 0-jet exclusive cross section and add events with more emissions below the cut. In case of a single shower emission we require also that the resulting Φ 1 point is projectable onto Φ 0 , as these are the only configurations at this order which are included in eq. (46) (see also eq. (10)). Both these conditions can be implemented with a careful choice of the starting scale of the shower. The preservation of the cross section below the cut is then ensured by the unitarity of the shower evolution. In practice, we allow for a tiny spillover up to 5% above T cut 0 in order to smoothen the transition. Shower events originating from Φ 1 and Φ 2 events require instead more care. The reason is that we must now preserve the NNLL +NNLO accuracy of the T 0 spectrum. This means that the Φ 2 points produced after the first emission must be projectable onto Φ 1 using the T 0 -preserving map mentioned earlier and discussed in detail in Ref. [23]. These constraints are most simply implemented by performing the first emissions in Geneva (using the analytic form of the NLL Sudakov factor and phase space maps) and only thereafter letting the shower act as usual, subject to the single restriction T 2 (Φ N ) ≤ T 1 (Φ 2 ). Since it can be shown that the shower acting on the resulting Φ 2 events alters the accuracy of the T 0 distribution only beyond NNLL [23], in practice we apply this procedure only to the Φ 1 events. We find that ) .
By choosing Λ 1 ∼ Λ QCD , the Sudakov factor U 1 (T cut 1 , Λ 1 ) becomes vanishingly small and we can relax the shower conditions on the 1-jet contributions. The showered events therefore originate from either dσ mc 0 or dσ mc ≥2 . The main differences which occur with respect to the Drell-Yan case are twofold: the decay of the Higgs boson and the choice of the starting scale for the gluon-fusion contributions. In our study, we have hitherto considered a stable Higgs boson. As the Higgs boson is a scalar with a small width over mass ratio, it is a legitimate strategy to separate the production process from the decay in the narrow width approximation. As long as we limit ourselves to a leading-order description of the decay process, it could be handled entirely by the parton shower program. In the following, we have chosen to work with a stable Higgs boson, in order to simplify the analysis routines, but we might as well have let Pythia8 decay it in order to achieve a more realistic description of the final state. The gluon-fusion contribution could be considered as giving rise to Φ gg 0 events which are different in nature to the standard Φ 0 events entering the Geneva formula. As explained in sec. 2 B 2, the reason is that the gluon-fusion contributions are entirely nonsingular and are therefore merely added at FO, lacking any resummed contribution. Since they enter at O(α 2 s ), any emission added by the shower would contribute beyond the claimed NNLO accuracy. We have, therefore, a greater freedom in the choice of the shower starting scale, which we set at the kinematic limit determined by the available centre-of-mass energy.

VALIDATION
In order to validate the FO accuracy of our results we have compared with a custom version of the Matrix code [65] which implements the Higgsstrahlung process. This program calculates cross sections at NNLO accuracy in QCD through its fully general implementation of the q T -subtraction formalism [66][67][68], in combination with the dipole-subtraction formalism [58,59] to deal with NLO-like singularities. To eliminate the dependence on the slicing parameter r cut ≡ q cut T,V H /M V H we numerically approach the limit r cut → 0 by the extrapolation procedure of Ref. [65] not only for inclusive cross sections, but also on a bin-wise level for distributions, as introduced in Ref. [69]. The uncertainties associated with this extrapolation procedure are combined with the statistical uncertainties to provide the overall numerical error of the predictions that is shown in figs. 5-7 and 9.
Our calculations are performed for pp collisions at the 13 TeV LHC using the PDF4LHC15 nnlo 100 PDF sets [70] available via the Lhapdf interface [71]. In the case of ZH production we restrict the invari- In fig. 5 we show the transverse-momentum spectra of the produced vector boson for each of the three possible V H processes, excluding for the moment the gluoninitiated channel in ZH production. The red hashed band associated with the Matrix curve reflects the error obtained from a simultaneous variation of the renormalisation and factorisation scales around the central scale (3point variation), while the blue band associated with the Geneva prediction has been obtained by following the procedure detailed in sec. 2 B 3. We observe good agreement for the central value in each case. The scale variation band in the Geneva case becomes slightly larger in the hard region at large transverse momentum because the mechanism for recovering the exact NNLO cross sections and scale variations discussed in sec. 2 B 3 becomes less effective, being based on the total cross section and not on differential distributions. From here on we leave behind the W ± H process and focus only on the ZH case, which displays broadly similar behaviour but exhibits a few more interesting subtleties due to the presence of the aforementioned gg production channel at NNLO. In figs. 6-7 we show comparisons of transverse-momentum and rapidity distributions between Geneva and Matrix, again neglecting the gg channel for the time being. At the value of T cut 0 chosen, we observe very good agreement with the Matrix results with the exception of the case of the Higgs boson transverse momentum. We find that agreement for this distribution is improved at a lower value of the cut at the expense of a higher statistical uncertainty (as previously discussed in sec. 2 B 4).
We now consider the inclusion of the gg channel at FO. In fig. 8 we show the impact at the partonic level in Geneva, focusing on the Higgs boson transverse momentum and the invariant mass of the V H system. We observe an effect of up to ∼ 20% on the differential distributions, demonstrating the importance of including this channel. We see also an increase in the scale uncertainties, related to the fact that the process gg → ZH is included in effect only at leading order (albeit O(α 2 s )). In fig. 9 we compare the Geneva predictions including the gg channel with those of Matrix for the Higgs boson transverse momentum and rapidity and again find good agreement.
Unfortunately, we were unable to compare with the results of similar calculations presented in Refs. [25,26]. The reason for this was that in the case of W ± H production, the authors in Ref. [25] neglected all top-quark effects while using a 3 × 3 CKM matrix. While each of these options can be set separately in the current public release of OpenLoops [73] which we rely on to provide amplitudes, the combination of both is not possible. In the case of ZH production, the same authors in Refs. [26] included the decay of the Higgs boson into a bb pair at NLO, which at present is omitted in our calculation. We intend to include the higher-order corrections to the decay in a future release of the program which will enable us to make a detailed comparison of the two results.

RESULTS
We now present our predictions for various spectra after interfacing with the parton shower provided by Pythia8 v8.235 [74,75]. For definiteness, we have chosen Pythia8's tune 18, we have set p ref T0 = 2.4 GeV, and have run with all matrix element corrections switched off, since now the radiative effects entering at higher order are provided by Geneva. In order to keep the analysis as simple as possible, we have also switched off all QED effects in the showering. In the following, we adopt the scale choice µ FO = M T V H but otherwise use the same values of the parameters as in sec. 3. Our scale uncertainty bands are calculated using a different procedure depending on whether quantities are either exclusive or inclusive in the additional radiation, as described in sec. 2 B 3. We reconstruct jets using the FastJet algorithm [76,77] with a jet radius R = 0.4 and a minimum p j,cut T = 30 GeV.
In fig. 10 we show the T 0 distribution at the partonic and showered level for ZH production with the gg channel switched off in three different regions: the peak (left pane), where the resummation effects are expected to be dominant; the transition region (centre pane) where the resummed and FO contributions should be on the same footing; and the tail (right pane), where the resummation is switched off and the FO perturbative expansion is a valid approximation. We confirm that, as expected, the NNLL accuracy of the T 0 distribution is preserved by the shower above T cut 0 . The shape below T cut 0 is determined entirely by Pythia8, but the cross section falling below the cut is preserved as required (apart from the small spillover discussed in sec. 2 B 6). The small contri- bution appearing between 0 < T 0 < T cut 0 at the partonic level is due to the nonsingular Φ 1 events which cannot be projected on a valid Born-like configuration and are therefore included only at fixed order.
In fig. 11 we turn on the hadronisation and show its impact on the showered distribution: we observe a large difference only in the peak region, as expected, with the corrections at larger values of T 0 being suppressed as O(Λ QCD /Q).
We continue with an examination of the effects of the shower on distributions other than T 0 . In fig. 12 we consider the transverse momentum of the Higgs boson and the rapidity of the hardest lepton, both quantities which are inclusive over any additional radiation. We note that showering does not significantly change the normalisation or the shape of the distributions, demonstrating that the NNLO accuracy is maintained even after showering and hadronisation. Additionally, we note that the scale variation uncertainties are unaffected by the shower and hadronisation stages. This is to be expected and a consequence of the fact that in this simple analysis we have neglected any uncertainty originating from the interface of our partonic predictions to the shower and from the hadronisation model. One could explore these effects in more detail by studying, for example, the variation in our predictions after modification of the shower starting scale for the former or explore different tuning parameters for the latter. Such investigation is beyond the scope of the current study.
In fig. 13 we show instead quantities exclusive in the additional radiation. Although we cannot claim NNLL accuracy in the resummation for these observables, we may anticipate that a certain amount of the accuracy from the prediction of T 0 may be inherited by other quantities. In the case of the rapidity of the hardest jet we see that the shower causes an overall shift of the distribution downwards by O(10%), most likely due to the acceptance cut on the jet above p j,cut T . Considering instead the transverse momentum of the V H system, we see that the shape of the distribution in the resummation region is significantly modified by the shower.
We now proceed to study the effect of the shower after including the gluon-fusion channel. The majority of distributions show similar effects, with the inclusive quantities surviving the shower stage unmodified. Two dis-tributions for which more significant differences are seen, however, are shown in fig. 14, where we plot the transverse momenta of the Higgs boson and of the V H system. We notice a significant deviation after the inclusion of the shower in the hard region. This is most likely to be a consequence of our choice of starting the shower at a very high scale for these contributions, and is in accordance with previous observations for similar gg-initiated processes [78,79]. Since the showering of these contributions starts at O(α 3 s ), in the present calculation we lack any further means by which we may constrain its effects. This motivates a future inclusion of the gg-initiated process at NLO.
Finally, in fig. 15 we examine the difference in the spec-tra when multi-parton interaction (MPI) effects are included in predictions for ZH production at the hadronised level and with the gg channel switched on. We observe that the T 0 distribution is significantly modified by inclusion of MPI, as already seen in Ref.
[30] -this follows from the definition of the beam thrust at the analysis level, which involves a sum over all final particles including those arising from secondary collisions. In the case of an inclusive quantity however (for example the Higgs boson transverse momentum), the shape of the distribution changes very little.

CONCLUSIONS
We have implemented the Higgsstrahlung process in the Geneva framework, which provides resummed predictions matched to the fixed-order calculation and a parton shower at NNLL +NNLO accuracy. In order to make a consistent choice for the profile functions used for the determination of the resummation scales, we have studied the interplay between the singular and nonsingular contributions in different regions of the Born-like phase space. As expected, we find that the region in which the resummation is applicable depends mostly on the value of the invariant mass of the V H system and it can be defined in a similar way for all classes of events showing only a mild dependence on the other kinematical variables.
We have confirmed the fixed-order accuracy of our results by comparison with the program Matrix and found very good agreement at T cut 0 = 1 GeV. Only for the Higgs boson transverse-momentum distribution do we find that an improved agreement can be reached by lowering the 0-jet resolution cutoff to T cut 0 = 0.1 GeV, at the price of an increased statistical uncertainty. This leads us to conclude that the discrepancy is most likely due to the missing O(α 2 s ) power corrections in Geneva as compared to Matrix.
We have provided predictions at the showered and hadronised levels by interfacing with the parton shower program Pythia8. We first confirmed that the accuracy of the T 0 distribution is unaffected by the showering and then showed that the inclusive distributions retain their NNLO accuracy and are mostly unchanged by the shower. The shower effects were found to be more significant for more exclusive distributions.
We have also included the gluon-fusion channel and studied the differential distributions which are affected most by its inclusion. We observe larger shower effects connected with these configurations, possibly related to the higher starting scale used for the shower.
Finally we were also able to include MPI effects and showed their impact on inclusive and MPI-sensitive distributions.
The code used for this study is available upon request to the authors and will be made public in a future Geneva release at http://geneva.physics.lbl.gov.
There are some clear directions in which this work could be furthered. At present, the decay of the Higgs boson in our implementation can only be provided at LO by Pythia8. A full NNLO calculation of the Higgs boson decay to a bb pair would be desirable in order to improve the description of the final states which are experimentally accessible. The combination of the production and decay processes in the narrow-width approximation is in principle feasible within the Geneva framework, and we plan to study this in a future work. Another avenue worth pursuing is the inclusion of subleading power corrections at the fully differential level, which would reduce the size of the neglected terms and thus improve the predictions for distributions even when a larger value of the resolution cutoff is used. These effects are likely to become more important as processes with more complex final-state phase spaces are considered.
Given the difficulty of discovering New Physics at the LHC, it is now more important than ever to be able to make precise predictions of the SM backgrounds both at the fiducial cross section level and when extrapolated over the full phase space. Since Monte Carlo event generators are the primary tool used to provide these predictions, it is vital that they are made as accurate as possible. This allows state-of-the-art theoretical calculations to be made available to experimental collaborations so that they can be used directly in analyses.