$VH+\text{jet}$ production in hadron-hadron collisions up to order $\alpha_s^3$ in perturbative QCD

We present precise predictions for the hadronic production of an on-shell Higgs boson in association with a leptonically decaying gauge boson and a jet up to order $\alpha_s^3$. We include the complete set of NNLO QCD corrections to both charged- and neutral-current Drell-Yan type contributions, as well as the previously known leading heavy quark loop induced contributions which involve a direct Higgs-quark coupling. As an application, we study a range of differential observables in proton-proton collisions at $\sqrt{s} = 13~\text{TeV}$ for both the charged- and neutral-current production modes. For each Higgs production process, we assess the improvement in the theoretical uncertainty for both the exclusive ($n_{\text{jet}} = 1$) and inclusive ($n_{\text{jet}} \geq 1$) jet categories. We find that the inclusion of the NNLO corrections to the Drell-Yan type contributions is essential in stabilising the predictions and in reducing the theoretical uncertainty for both inclusive and exclusive jet production for all three modes. This is particularly true in the kinematical regimes associated with low to medium values of the transverse momentum of the produced vector boson and where the differential cross sections are the largest. For the neutral-current process, we find that the heavy quark loop induced contributions have their largest phenomenological impact (an increase in the size of the NNLO corrections, a distortion of the distribution shape and an enlargement of the left over remaining uncertainties) in kinematical regions associated to large values of $p_{T,Z}$ (typically above $150~\text{GeV}$) where the cross sections are smaller.

Since the discovery of the Higgs boson [1,2], the continued measurement of the properties and interactions of the Higgs boson with other fundamental particles has been a main goal of the LHC physics programme. In particular, this concerns the coupling strengths of the Higgs boson to known Standard Model (SM) particles such as gauge bosons and fermions through the study of differential observables related to various Higgs boson production and decays channels.
A key channel is the associated production of a Higgs boson with an electroweak gauge boson, referred to as the Higgs-Strahlung or V H channel (V = W ± or Z). The presence of a leptonically decaying gauge boson in the final state offers a handle to efficiently reject the contribution of multi-jet backgrounds. This is critical for the measurement of the hadronic decays of the Higgs boson, which are otherwise overwhelmed by such background processes. An important example is the decay of the Higgs boson into a pair of b-quarks, which in the SM occurs with an expected branching fraction of ≈ 58%. Measurements of this decay channel (which are uniquely accessible in the V H channel) are essential to directly probe the interaction strength of the Higgs boson with quarks, and also provide important constraints on the total decay width of the Higgs boson.
Experimentally, the H → bb decay was observed by the ATLAS [3] and CMS [4] collaborations in 2017, precisely through the V H production mode with a significance of 5.6 and 5.3 standard deviations for CMS and ATLAS respectively. In addition, first differential measurements based on simplified template cross sections as a function of the transverse momentum of the vector boson have been reported in ref. [5] and updated including the full Run II data set in refs. [6,7]. The latter measurements indicate that the observed production rates are presently consistent with SM expectations within experimental uncertainties of the order of 20 percent. These measurements are performed by categorising the observed event rates according to the number of charged leptons ( ) and hadronic jets. Such selections are necessary to improve the experimental sensitivity to the signal processes. In all channels, events are required to have exactly two b-tagged jets, which form the Higgs boson candidate.
The 0, 1 and 2 lepton channels respectively provide access to the V H decay processes: ZH → ννbb, WH → ν bb and ZH → + − bb. Events are then further categorised based on the presence of hadronic jets which are not b-tagged (i.e. in addition to those two which form the Higgs candidate). Examples are: n jet = 0, no additional hadronic jets; n jet = 1, exactly one additional hadronic jet; n jet ≥ 1, at least one additional hadronic jet. Notably, almost half of the selected events in an experimental analysis [6,7] are contained in the n jet = 1 or n jet ≥ 1 categories. These latter categories will be the focus of this work.
As the experimental analyses rely on the theoretical knowledge of the V H (+jet) production and decay modes to interpret the data, it is of crucial importance to have the most accurate theoretical predictions available for differential observables in the selected event categories. Depending on the event category, various levels of theoretical precision are currently available: At fixed-order accuracy, NNLO QCD corrections for the V H production mode have been available for some time [8][9][10][11][12]. More recently, calculations which also include decay corrections of the Higgs boson up to NLO [13] and NNLO QCD accuracy in the five-flavour [14][15][16] 1 and four-flavour schemes [18,19] have been presented. Other NNLO QCD calculations including an interface to a parton shower event generator have been presented in [20,21] and [22], where the latter includes decay corrections. The phenomenological impact of heavy quark loop contributions has been studied in [23][24][25][26][27]. For the V H + jet production mode, fixed-order NLO QCD computations for V H and V H + jet production have been merged in the context of parton showers to provide full NLO+PS simulations for these (on-shell) Higgs production modes [22,28]. The corresponding results have been further complemented with electroweak corrections in [29]. More recently, predictions including NNLO QCD corrections for differential observables related to the W + H + jet production mode were presented (by us) in ref. [30].
It is the purpose of this paper to provide precise predictions for observables related to all three V H + jet production modes with V = Z, W + , W − where the vector boson decays leptonically and the Higgs boson is produced on its mass shell. In particular, we increase the perturbative accuracy of the Drell-Yan-type component of the cross-section (which is numerically dominant) to NNLO QCD accuracy for all production modes. The extension of these predictions to higher orders, including corrections to Drell-Yan type and heavy-quark loop-induced contributions up to the third order in α s , provides a new level of theoretical precision for these V H + jet production modes. This improvement will be vital for future interpretations of experimental measurements, as the current knowledge of the theoretical modelling of the signal process (currently limited to NLO accuracy in n jet > 0 categories) constitutes one of the dominant systematic uncertainties in signal extractions of these Higgs-Strahlung processes [31].
The structure of the paper is as follows: In section 2 we summarise the general structure of the calculation, and provide details of the various ingredients which enter the computations for the three distinct V H + jet production modes, up to order α 3 s . Our numerical results are presented in sections 4 and 5 for the case of pp collisions at 13 TeV, where both inclusive and exclusive jet categories are considered and compared. Specifically, we provide predictions for fiducial cross sections in section 4 and differential distributions in section 5, and discuss the phenomenological importance of including higher-order corrections in both cases. Section 6 presents a summary of the results obtained together with an outlook for further studies. Some of the results associated to W + H + jet production have been previously presented in ref. [30] and will be included here to facilitate a direct comparison with the other production modes, W − H + jet and ZH + jet, which are presented here for the first time.

Details of the calculation
In this work, differential observables for the production of a Higgs boson in association with a weak gauge boson and at least one resolved jet are computed up to and including order α 3 s contributions in perturbative QCD. More specifically, we consider the two charged-current processes pp → HW ± + jet → H + ± ν + jet (2.1a) and the neutral-current process, where the Higgs boson is produced on its mass shell and the weak boson decays leptonically.
The calculation is carried out within the NNLOJET framework [32] which is a fixed-order parton-level event generator that employs the antenna subtraction formalism [33][34][35][36][37][38][39][40][41] for the treatment of infrared divergences in order to retain the full kinematic information on the final state.
At leading-order, the processes in Eq. (2.1) give rise to cross-section contributions that are proportional to α s . These are given by Higgs-Strahlung diagrams, where a Higgs boson is emitted from a massive gauge boson of an underlying Drell-Yan type V + jet process.
An example of Born-level diagram is depicted in figure 1a. In the present calculation, full NNLO QCD corrections to the Drell-Yan type contributions (denoted as DY) are considered, i.e. up to order α 3 s for all three Higgs production modes. Starting from order α 2 s , an additional class of subprocesses contributes in which the Higgs boson couples to a virtual heavy quark loop. In our calculation, we consider leadingorder contributions from such heavy quark loop induced contributions up to order α 3 s : For both charged-and neutral-current processes this includes the top-quark loop corrections at order α 2 s y t , denoted as the R I -type contributions and depicted in figure 1b. In the case  at order α s for the Drell-Yan process, denoted as DY-type; (b) at order α 2 s for the top-loop induced processes, denoted as R I -type; and (c) at order α 3 s for the one-loop gluon-gluon initiated process gg → HZg, denoted as ggF-type. Details on the individual contributions given in the main text.
of the neutral-current process, additional loop-induced contributions related to the gluonfusion process gg → HZg must be considered. Those are proportional to α 3 s and illustrated in figure 1c. In the following those are denoted as ggF-type.
The various ingredients necessary for the computation of the W + H + jet mode have been briefly outlined in ref. [30]. Similar ingredients are included in the computation of observables associated to the W − H + jet and ZH + jet production modes and presented here for the first time.
In the following, the relevant details for the calculation of the DY-type and heavy quark loop induced contributions are provided in sections 2.1 and 2.2 respectively while section 2.3 presents a summary of the cross section contributions for the different production modes, up to α 3 s .

Drell-Yan-type contributions
The Drell-Yan-type contributions are characterized by a Higgs boson coupling to an intermediate weak gauge boson and are thus proportional to g 2 V V H at the cross-section level. They are denoted as dσ V, DY in the following with V ∈ {W ± , Z} and are computed up to order α 3 s in this work. As evident from Eq. (2.1), we consider leptonic final states from the off-shell weak gauge boson. The amplitudes for this set of contributions can be obtained from the corresponding Drell-Yan process without the Higgs boson, provided that external momentum conservation has not been assumed in their derivation. In that case, the off-shell W ± or Z boson can be substituted by a gauge boson that subsequently emits an on-shell Higgs particle and then decays into leptons. The necessary Drell-Yan amplitudes for all three V H + jet modes up to order α 3 s can therefore be constructed from the various components of the existing calculations for the W + jet [42] and Z + jet [32] processes within NNLOJET.
The above construction was directly applicable in most cases. In a few cases, new analytic results for the amplitudes that did not impose momentum conservation on the external states were derived first. In particular, the tree-level and two-loop amplitudes, which constitute the building blocks at the double-real and double-virtual levels to the partonic cross section were obtained in this manner. At the real-virtual level, all one-loop amplitudes are obtained through the OpenLoops 2 [43,44] library (see also [45]). This library was further used to validate the manually constructed tree-level amplitudes as well as to evaluate the four-quark tree-level amplitudes in the final computation.
As the infrared structure of the squared amplitudes is independent of the gauge-boson species, the construction of subtraction terms to treat infrared singularities is completely analogous between the various production modes of the Drell-Yan category. Furthermore, after appropriate substitution of the colour-ordered matrix elements, all subtraction terms for the computation of the Drell-Yan-type contributions to the V H+jet process up to order α 3 s can be constructed from the existing V + jet calculations [32,42] (see references therein for the appropriate antenna functions).

Heavy quark loop induced contributions
Starting from order α 2 s , an additional class of corrections contribute that are characterised by a Higgs boson coupling to a heavy quark loop. These can be treated independently from the Drell-Yan-type discussed above, as the corresponding amplitudes are separately gauge invariant. We have included the leading-order heavy quark loop induced contributions up to α 3 s for all three production modes, using the process libraries for the relevant squared amplitudes as provided by OpenLoops 2 [43]. These contributions are both ultraviolet and infrared finite and therefore do not require any dedicated subtraction procedure.
At order α 2 s , both the charged-current and neutral-current processes receive a contribution that is characterised by the Higgs boson (but not the gauge boson) coupling to the closed quark loop-commonly referred to as R I -type amplitudes [23]. The corresponding cross section contribution is given by the interference of the tree-level Drell-Yan-type amplitude with the R I -type amplitude shown respectively in figures 1(a,b). To this end, we consider the numerically dominant top-quark loop contribution that is proportional to α 2 s y t and denote its contribution to the V H + jet cross section by dσ V, R I . From the phenomenological point of view, taking into account only the leading-order R I -contribution is found to be sufficient given its relatively small numerical impact, as also pointed out in ref. [30] for the case of W + H + jet production. 2 In the case of the neutral-current process, additional loop-induced contributions are included that belong to the gluon-fusion channel gg → HZg and which are proportional to While formally suppressed by a factor of α s with respect to the R I contributions, the large gluon luminosity in LHC proton-proton collisions strongly enhances the ggF-type corrections such that its phenomenological impact can be even more sizeable than the R Itype terms in ZH + jet production. They further lead to distinct features in the fiducial cross section and differential distributions (as compared to the charged-current processes, where such contributions do not exist). These effects will be elaborated on in sections 4 2 The extension of this calculation to higher order relies on accounting for all gauge-invariant sets of fermion loop contributions at a given order in αs. The computation of the order α 3 s contribution relies on the knowledge of a set of heavy quark-loop amplitudes which are only partially known. An illustrative example of such unknown contributions for a quark-antiquark initiated reaction is: t and 5.

Cross section contributions up to order α 3 s
We conclude this section by summarizing the specific cross section contributions entering our final (N)NLO predictions. To this end, we use a labelling based on the underlying Drell-Yan-type contributions, for which we follow the customary perturbative counting in terms of the strong coupling: The R I -type contributions from heavy quark loops that start to contribute from O α 2 s are consistently included (at leading order) together with NLO corrections to the DY-type contributions. While the ggF-type contributions are formally suppressed by an additional power of α s compared to the R I -type contributions, we choose to include them also starting from NLO. The reason for this is twofold: The large gluon luminosity at the LHC enhances this contribution to the extent that it is of similar size to both the NLO DY-type and R Itype corrections. Furthermore, the combination of the R I -and ggF-type contributions together with the NLO DY-type corrections corresponds to the state of the art prior to our work [22,26,28,29] and thus facilitates a direct comparison of the numerical impact of the newly computed NNLO corrections to the DY-type parts.
The labelling of (N)NLO predictions follows the definitions below: and will be referred to throughout the remainder of this paper to specify the contributions included in our numerical results for fiducial cross sections and kinematic distributions.
In addition, we also provide predictions that maintain a strict power counting in α s when presenting the fiducial cross section results in Tables 1-3. To this end, we define the , which represents the NLO cross section that only includes terms up to O α 2 s , i.e. without the ggF-type contribution for the neutral-current process. For consistency, we use this notation in all three tables although in the charged- In the neutral-current production mode, this enables us to quantify the impact of the inclusion of the α 3 s corrections arising either from the ggFtype contributions or from the pure DY-type contributions, as will be further elaborated on in section 4.
-8 -In this section we review the calculational set-up as well as the kinematical constraints that are imposed to define the fiducial cross sections for the class of V H + jet processes. We consider the associated production for all possible massive gauge bosons, V ∈ {W ± , Z}, and separately report the 1-jet inclusive (n jet ≥ 1) and 1-jet exclusive (n jet = 1) cases.
Throughout the paper we will refer to these latter selections as inclusive and exclusive respectively, which strictly denote the jet multiplicity selection.
We provide predictions for proton-proton collisions at √ s = 13 TeV using the parton distribution function NNPDF31_nnlo_as_0118 [46] from the LHAPDF library [47]. The inclusive (exclusive) processes are required to contain at least (exactly) one jet with transverse momentum p ⊥ > 20 GeV. These requirements define the inclusive and exclusive production cross sections, which are respectively denoted as σ ≥1j and σ 1j . The anti-k t jet algorithm with the parameter ∆R = 0.5 is used to cluster final-state partons into jets.
Charged leptons are required to have a transverse momentum p ⊥, ± > 25 GeV with the rapidity constraint |y ± | < 2.5. In case of the charged-current processes, W ± H + jet, an additional requirement on the minimum missing transverse energy larger than 25 GeV is imposed for the neutrinos.
For the electroweak input parameters, we employ the G µ -scheme with the full list of independent parameters entering the computation given by The running of the strong coupling (α s ) is evaluated using the LHAPDF library with the associated PDF set. Finally, in the case of W ± H production, a diagonal CKM matrix is used for the vector-boson-quark couplings. We make this approximation based on the expectation that the phenomenological impact from non-diagonal effects are sub-leading for the processes under consideration. In order to estimate theoretical uncertainties from missing higher orders, we follow the conventional 7-point variation prescription: with the constraint 1 2 ≤ µ F /µ R ≤ 2. All inclusive predictions adopt this error estimate for theory uncertainties.
In the case of exclusive predictions, the cross section can also be viewed as the difference between the one-and two-jet inclusive calculation (3.2) and the above prescription corresponds to a correlated variation of the scales in the individual inclusive predictions.
In the context of exclusive results, we will therefore often refer to the standard 7-point variation as the correlated error prescription.
As a more conservative estimate, we additionally consider a second method when presenting differential distributions that is based on ref. [48] and referred to as the uncorrelated prescription. Here, the exclusive cross section and uncertainty measure are given as: where ∆ ≥1 (2)

Fiducial cross sections
The cross-section predictions utilizing the fiducial cuts described in section 3 are summarised in tables 1-3 for W + H+jet, W − H+jet, and ZH+jet production, respectively. In these tables, the rows correspond to results for the inclusive (≥ 1 jet) (top row) and exclusive (1 jet  results where all contributions up to order α 3 s are taken into account and are denoted by σ NNLO .
Only the error estimates based on the 7-point scale variation, i.e. the correlated prescription described in section 3 are quoted for the sake of clarity. We now discuss the impact of including the NNLO DY-type corrections for the fiducial cross section for all three associated Higgs boson production mode in the inclusive and exclusive jet cases separately.

Inclusive predictions
We first focus our attention on the inclusive results presented in the top row of tables 1-3.
As expected, the fiducial cross section values which involve purely Drell-Yan like contributions, display well converging perturbative behaviour as higher and higher orders are considered. Indeed, by considering the numbers displayed in columns 2 and 4, we observe a notable stabilization of the inclusive cross section at NNLO: the correction to the cent-  The impact of the heavy quark loop contributions on the inclusive cross section can be studied by contrasting the results presented in the second and third columns (at NLO) and the fourth and fifth columns (at NNLO). Here, the results for charged and neutral channels start to differ significantly and we observe the following behaviour: WH + jet: The numerical impact of the top-loop amplitudes appearing at order α 2 s (NLO) is relatively small and of comparable magnitude to the Drell-Yan-like α 3 s (NNLO) contribution. This behaviour has already been noticed for the W + H + jet production mode studied in [30] and highlights that the inclusion of the heavy quark loop contribution is mandatory for precision phenomenology.
From this comparison, we can further infer that the α 2 s top-loop contributions are likely to be as robust to scale variations as the NNLO Drell-Yan-like predictions and that further higher-order corrections to this contribution are only relevant in conjunction with the Drell-Yan-type corrections at N 3 LO. The latter are however unknown and presently beyond the reach of perturbative computations.
ZH + jet: The heavy quark loop induced contributions in ZH + jet production enter not only at order α 2 s through the R I -type corrections, similarly to the charged-current case, but also at order α 3 s through the gluon-gluon-initiated ggF-contributions. The impact of heavy quark loop corrections is thus noticeably different for this Higgs boson production mode at the different orders in α s . At O α 2 s , the R I -type terms contribute to a slight increase of the cross section but have no sizeable effect on the uncertainties. The impact of the gluon-gluon-induced α 3 s terms is positive, large in magnitude and contributes to the cross section with a residual uncertainty that is enlarged compared to the case when only Drell-Yan like contributions are included.
As noted in section 2, the impact of the gluon-gluon-induced contributions is significantly enhanced due to the high gluon luminosity of 13 TeV proton-proton collisions.
These contributions form a gauge-invariant subset of the α 3 s corrections and their inclusion is justified and in fact necessary for the correct phenomenology of the process.
Numerically, they provide the bulk of the order α 3 s corrections (i.e. they dominate over the Drell-Yan type corrections at this order), which overall leads to an O(17%) increase compared to the previous order and a much increased theoretical uncertainty of O(5%).

Exclusive predictions
Next, we consider the exclusive fiducial cross sections given in the bottom row of the We reiterate that the results shown in tables 1-3 show uncertainties obtained with the correlated scale prescription. This is noted here because there is a partial cancellation in uncertainties for the exclusive ZH + jet case at the NNLO level, i.e. the last entry in the second row of table 3. If instead the correlated scale prescription is used, the uncertainty for ZH + jet at NNLO is ∆ NNLO uncorrelated = 0.42 fb. The impact of including the heavy quark loop induced contributions at the differential level, in terms of size/shape of the corrections and residual scale uncertainty at NNLO (including scale uncertainty prescriptions) for all three Higgs boson production modes, will be investigated for both inclusive and exclusive selections in section 5.
-13 -In this section we turn our attention to differential observables for the three V H + jet processes as described in section 2 and using the numerical set-up given in section 3 for the 13 TeV LHC. We focus in particular on observables that are relevant for the determination of the V H signal in experimental analyses [6]. Accordingly, we choose to separate the presentation of our results in three parts. First, in section 5.1, we compare the inclusive and exclusive Higgs boson production modes for the transverse momentum spectra of Higgs boson and vector boson for the W − H + jet, and ZH + jet processes. Second, and motivated by the necessity to suppress overwhelming multi-jet backgrounds from tt production where one of the top quarks decays semi-leptonically, we present results for the exclusive charged current production modes in section 5.2. Third, we discuss results for the inclusive neutral current production mode in section 5.3 because typically any number of jets is accepted in the experimental analysis to increase the signal acceptance.
For all results presented in section 5, the computation of the differential cross sections is performed using the expressions presented in section 2.3. More specifically, the (N)NLO expressions are given in eq. (2.3a) for the charged current case and in eq. (2.3b) for the neutral current process. In other words, the previously known heavy quark R I -type contributions starting at order α 2 s and, in case of the neutral-current process, the gluon fusion contribution appearing first at O(α 3 s ) are included in both NLO and NNLO predictions. With this choice, the study of the impact of the pure NNLO (i.e O(α 3 s )) DY-type contribution computed for the first time here is facilitated.

Inclusive vs. exclusive predictions
Our aim here is to compare the results obtained for predictions in the 1-jet inclusive (n jet ≥ 1) and 1-jet exclusive (n jet = 1) categories at a given fixed order in α s . We shall here focus on the transverse momentum spectra of the Higgs boson and the vector boson in the production processes W − H + jet, and ZH + jet and highlight similarities and differences for these two production modes. We begin our discussion with a general description of the figures illustrating these differential results. Specifically, in section 5.      towards higher values of p T , indicating that the high transverse momentum region is dominated by VH +2-jet production, which is present at NLO level only in these computations.
In both exclusive and inclusive cases, the inclusion of the O(α 3 s ) DY-type corrections, stabilises the predictions for the transverse momenta spectra, and leads to a reduction of the residual theoretical uncertainties. This is most pronounced in the exclusive production mode, where the inclusion of the NNLO corrections has a bigger impact.

ZH + jet distributions
For the distributions for the neutral production mode presented in figure 3, we observe that the impact of the inclusion of the NNLO corrections has a similar qualitative effect on the inclusive and exclusive K-factors. As in the charged-current case, we observe that the size of these K-factors in the exclusive category is larger than in the inclusive category. This is clearly seen in the second and third panels, together with the veto efficiency in the fourth panel.
In addition, we notice a distortion of the shape of the K-factors for both distributions which is absent in the charged current mode. This is present in both inclusive and exclusive The impact of the DY-and non-DY-type corrections present at O(α 3 s ) will be further analysed in section 5.3 for a wider set of observables in the inclusive neutral production mode.

The exclusive charged-current processes
In this section, we focus on the 1-jet exclusive WH process. As previously noted, the motivation for considering a 1-jet exclusive selection is that the large background from the semi-leptonictt channel can be suppressed to a large degree. Our aim is to compare and contrast the effect of the NNLO corrections on the differential predictions for the two charged-current processes. We consider four observables: the transverse momenta spectra of the vector boson and the Higgs boson, as well as the transverse momentum and rapidity distribution of the jet shown in figures 4-7. These figures contain two panels: The top panel shows the absolute distributions at each perturbative order in α s and the lower panel displays the respective K-factors. In the latter, the denominator of the NLO ratio contains DY type and heavy quark loop R I type contributions up to order α 2 s as detailed in section 2.3.

The transverse momentum spectra of vector and Higgs bosons
The behaviour of both the gauge boson and the Higgs boson transverse momentum spectra for the W − H + jet process were analysed in section 5.  To compare the two production modes, we therefore focus on the K-factors presented in the second panels of figures 4 and 5. We see that the corrections are sizeable and only the more conservative uncorrelated error estimate is able to provide overlapping uncertainty bands when going from NLO to NNLO. We further observe that the NNLO corrections and the associated uncertainties at NNLO are always larger in the W − case (up to 20%) compared to the W + case (which are about 10%). A possible interpretation is that the relative contribution of the gg and qq partonic channels is larger for the W − mode, and this results in larger theoretical uncertainties. This is particularly true for the regions of kinematic distributions which are sensitive to the input PDFs at large-x values (i.e. forward rapidities), where the d/u ratio becomes small and the impact of these NLO channels increases.

Jet related observables
We now consider predictions for two specific jet related exclusive observables, namely, the transverse momentum and the rapidity of the leading jet presented in figures 6 and 7 respectively.
For both of these exclusive Higgs boson production modes, the inclusion of the NNLO corrections lead to a stabilisation of the predictions with reduced left-over theoretical uncertainty bands compared to the results obtained at NLO. However, the NNLO corrections The results for the rapidity distribution of the leading jet shown in figure 7 also exhibit some marked differences between the W + and W − production modes. The most noticeable differences occur in the forward rapidity region, where larger NNLO K-factors and residual uncertainties in the W − case are observed.
In both cases, the NNLO K-factors relative to NLO, are negative for centrally produced jets, −12% for W + and −18% for W − . However, the NNLO corrections are largely independent of y j 1 for the positive charged-current process, while for the W − H+jet process, the NNLO corrections are small at large |y j 1 |. In both cases, the residual scale uncertainties are considerably reduced at NNLO, and are typically ±6% for the production of a W + and ±8% for the W − case.
The difference in behaviour can be attributed to the parton distribution functions and the differences in the luminosity of the partons that initiate the hard scattering. At large rapidities, high momentum fractions x are probed that can vary significantly between the two charged-current processes where valence-and sea-quark contributions are interchanged.
Note that similar qualitative features were also observed when comparing the W + + jet and W − + jet processes (with no Higgs boson in the final state) for the leading-jet rapidity distributions in ref. [42].

The inclusive neutral-current process
In this section, we present differential predictions related to the associated production of a Higgs boson, a neutral Z boson and at least one hadronic jet in the final state, i.e 1-jet inclusive production (n jet ≥ 1). We focus here on three categories of inclusive

Jet observables in ZH +jet production
Focussing first on the leading jet p T spectrum displayed in figure 9a, we observe that the NNLO corrections are (a) mostly independent of p T,j and at the percent level, (b) clearly lead to a stabilization of the perturbative predictions with scale uncertainties that are contained within the NLO uncertainty band, and (c) have residual scale uncertainties of about ±8% compared to NLO, over the whole kinematical range.
The leading jet rapidity distribution is shown in figure 9b. We see that the size of NNLO corrections and corresponding scale uncertainty band are considerably larger than those observed for the leading jet transverse momentum observable. From the second panel, we observe that compared to NLO, the NNLO prediction is small in the central rapidity region but increases with |y j 1 |, with little reduction of scale uncertainty band.
Comparing the NLO and NLO DY results leads us to identify the inclusion of ggF-type  loop contributions as significant (about +20% compared to NLO) and uniform in p T,j and y j 1 . The inclusion of the ggF contribution changes the phenomenological picture at NNLO specially in the most forward rapidity region, where the cross section is the smallest.

ZH +jet specific discussion
Finally, we discuss results that are specific to the neutral-current Higgs boson production process. To this end, we consider observables that are motivated by the experimental measurements and used in multivariate discriminants to improve the sensitivity of the ZH analyses [6]. Figure 10a shows the invariant mass distribution of the Higgs boson-jet system in ZH + jet inclusive production. This observable is important for the rejection of multi-jet backgrounds. For the production of a Higgs boson and a jet of p j T > 20 GeV and above the kinematic threshold of M Hj 1 ∼ 150 GeV, where a maximum is reached, the spectrum is steeply falling. Focussing on the second panel of this figure 10a, we observe that the relative size of the NNLO corrections increases from −10% at M Hj 1 ∼ 100 GeV to +30% at M Hj 1 ∼ 750 GeV with residual scale uncertainty bands increasing from ±3% to ±10% over the same range. This behaviour can be attributed to the fact that at high invariant mass, multi-jet configurations, included here at NLO level only, dominate compared to the exclusive 1-jet production process, which is here clearly suppressed.
The angular separation between the Higgs boson and the Z boson is an important input to identify the Higgs boson candidate and can be probed via the azimuthal angle and the pseudo-rapidity separation variables shown in figures 10b and 10c, respectively.
As argued in section 5.1, the Z and H bosons are preferably produced back-to-back in the transverse plane. This statement is further corroborated by the azimuthal angle distribution shown in figure 10b that strongly peaks around ∆φ ZH ∼ π. While additional emissions further open up the phase space for ∆φ ZH → 0, the region below this peak value still remains severely suppressed given the large mass of the two objects. For values of the azimuthal angle separation close to zero, the cross section value is also close to zero. This is clearly seen in the first panel of this figure 10b. To avoid unstable bin to bin fluctuations, for both of the panels in this particular figure 10b the size of the corrections in the first two bins has been averaged. Looking at the second panel, we see here that, unlike in all the distributions discussed so far in this inclusive neutral-current mode, the size of the corrections is the largest where the cross section is the largest. Indeed, close to the peak value of ∆φ ZH ∼ π, attained exactly in a leading order back-to-back configuration for ZH production, the size of the NNLO corrections relative to NLO and the corresponding uncertainties are at the level of 10% while they decrease to be at percent level for ∆φ ZH → 0. This observable is sensitive to the polarisation of the Z boson and helps in distinguishing the ZH signal from the Z + jet backgrounds. Being defined in the rest frame of the gauge boson and considering the fact that the leptonic decay Z → − + does not take part in any QCD interactions, the observable cos(θ CS ) is largely insensitive to the recoil from the jet and other QCD emissions. As a consequence, QCD corrections to this observable lead to remarkably flat K-factors (close to one), that are essentially given by those of the fiducial cross sections given in table 3. This is best seen comparing the results labelled as (N)NLO and NLO DY in the second panel of figure 10d.
For the majority of the distributions presented here, we observe that at low or medium values of the kinematical variable probed and corresponding to a region where the differential cross sections are the largest, the impact of the NNLO corrections are mainly due to genuinely new Drell-Yan type contributions present at order α 3 s . The inclusion of those corrections lead to a clear stabilisation of the predictions and a considerable reduction of the left over theoretical uncertainty band. The inclusion of the quark-loop induced contributions is relevant to achieve percent-level accuracy in those regions. Instead, at higher values of the kinematical variable probed, the impact of the quark-loop induced contributions becomes significantly more important and one observes a magnification of the size of the corrections, a distortion of the distribution shape, and an enlargement of the left-over remaining uncertainties in the NNLO predictions.

Conclusions and Outlook
In this work we have presented state-of-the-art predictions for the charged-and neutralcurrent V H + jet production processes. These include the corrections to the Drell-Yan type contributions up to order α 3 s , i.e. at NNLO in perturbative QCD, as well as the leading heavy quark loop induced contributions up to order α 3 s . These calculations were used to perform phenomenological studies of all three V H + jet production modes in protonproton collisions at √ s = 13 TeV. In particular, a range of differential observables have been considered both for exclusive (n jet = 1) and inclusive (n jet ≥ 1) jet production.
One of the motivations for considering V H + jet production is that event selection categories in the experimental analyses with n jet ≥ 1 can often receive the largest contribution of signal events in an experimental set-up [31]. Note that requiring a jet in association with the V H final state typically lowers the perturbative accuracy of a fixed order in α s prediction, which therefore increases the relative theoretical uncertainty associated with the modelling of the V H + jet signal in these categories (and overall reduces the sensitivity to interpret the Higgs data). It is for these reasons that dedicated calculations applied to the V H + jet categories (such as those presented here), are vital for future measurements/interpretations in the V H channel.
We have found that the inclusion of higher-order QCD corrections to the Drell-Yan type contributions are essential in stabilizing the predictions and in reducing the theoretical uncertainty for both charged-and neutral-current production modes (see also [30] where we presented selected results for W + H + jet production). This is particularly true in the kinematical regimes associated with low to medium values of the transverse momentum of the produced vector boson and where the differential cross sections are the largest.
The heavy quark loop induced contributions to the charged-current processes of order α 2 s y t introduce corrections at the percent level. While knowledge of the QCD correction to these contributions is desirable, the (known) leading contributions of order α 2 s y t are likely sufficient to obtain cross-section predictions accurate at the percent-level.
For the neutral-current process, we find that the heavy quark loop induced contributions have their largest phenomenological impact via the gluon-gluon fusion type contributions, named as ggF, and first present at order α 3 s . We find that, due to the inclusion of these quark-loop induced contributions in the computations one observes an increase in the size of the NNLO corrections, a distortion of the distribution shape, and an enlargement of the left over remaining uncertainties in the NNLO predictions. In most cases, however, these effects appear enhanced in kinematical regions associated to large values of p T,Z (typically above 150 GeV) where the cross sections are smaller. A reduction of this source of theoretical uncertainty will likely require the computation of order α 4 s corrections. The calculation of the necessary two-loop amplitudes enabling such a computation will be a challenging task.
Overall, the calculations presented in this work are critical for improving the precision of the Higgs boson signal templates [49,50], which play an essential role in the measurement and interpretation of Higgs boson production data at the LHC. In particular, these calculations are relevant for both the exclusive (n jet = 1) and inclusive (n jet ≥ 1) event selections, for both charged-and neutral-current processes in the V H production mode.
In the long term, the fixed-order calculations presented here also serve as an important milestone towards achieving an NNLO+PS accurate [51] description of the V H+jet process, deriving jet-vetoed cross sections at N 3 LO for associated Higgs production, and eventually aiming for a fully differential N 3 LO calculation for the V H process. Given the large data samples which are anticipated to be collected during the high-luminosity phase of the LHC, this level of theoretical precision is desirable to increase the sensitivity of future Higgs boson measurements.