text:
\section{Introduction}
\label{sec:intro}

After the discovery of the Higgs boson in
Run~I~\cite{Aad:2012tfa,Chatrchyan:2012xdj}, one of the main tasks of
the ongoing LHC~Run~II is to perform accurate measurements of Higgs
properties. In order to carry out this precision physics program, it is
important to study Higgs production in all the main production modes,
and compare measurements with theory predictions, for total
cross sections and differential distributions. An important goal which
is expected to be achievable with the Run~II full luminosity is to
establish solid statistical evidence for \textrm{HV}{} associated
production~\cite{Aaboud:2017xsd,Sirunyan:2017elk}.

The past years have seen a remarkable progress in NNLO QCD
calculations, and, currently, all $2\to 2$ SM scattering processes are
known to this accuracy, see~\emph{e.g.}
ref.~\cite{Heinrich:2017una}. Thanks to this progress, the description
of colour singlet final states has reached a high level of
accuracy. This is particularly true for processes where, at
leading order (LO), there are no gluons in the initial state: in this case
higher-order corrections are typically moderate, and hence including
NNLO corrections leads to very stable results, with small perturbative
uncertainties. For \textrm{HZ}{} production, the NNLO corrections have been
computed for the inclusive cross section~\cite{Brein:2003wg} as well
as for differential
distributions~\cite{Ferrera:2014lca,Campbell:2016jau,Caola:2017xuq}. Electroweak
corrections for this process are also known at NLO for inclusive
cross sections and differential
distributions~\cite{Ciccolini:2003jy,Denner:2011id}, and are
implemented in the public code {\tt
  HAWK}~\cite{Denner:2014cla}. Notably, NLO electroweak and QCD
corrections were simultaneously matched to a parton shower in
ref.~\cite{Granata:2017iod} for \textrm{HV}{} and \textrm{HV+jet}{}, using, in the
latter case, the \noun{MiNLO}{} method (to be described below).

Since associated \textrm{HV}{} production has a small cross section, it is
often considered in association with a Higgs boson decaying to a
$b$-quark pair, which is the largest Higgs decay mode. In this case
properties of the $b$-jets arising from the Higgs decay products are
used in experimental analysis to enhance the signal over SM
backgrounds, and they will also be important for extracting precise
information on the $b$-quark Yukawa coupling, especially in the Higgs
boosted regime. Because of this, including QCD corrections to the
${\textrm{H}}\!\to\!b\bar{b}{}$ decay is particularly important, especially since these
corrections are known to be large. The QCD NLO corrections to the
Higgs decay to massive $b$-quarks have been known for a long
time~\cite{Braaten:1980yq,Drees:1990dq,Sakai:1980fa,Janot:1989jf},
whereas more recently NNLO corrections were computed in
refs.~\cite{Anastasiou:2011qx,DelDuca:2015zqa} for massless
$b$-quarks. In the last few years the focus has moved towards a
combination of the aforementioned fully-differential NNLO computations
for $pp\to\textrm{HV}{}$ with differential NLO and NNLO results for
${\textrm{H}}\!\to\!b\bar{b}{}$. The current state-of-the-art results are those obtained in
refs.~\cite{Ferrera:2017zex,Caola:2017xuq}, where the
fully-differential QCD NNLO computations for $pp\to\textrm{HV}{}$ and ${\textrm{H}}\!\to\!b\bar{b}{}$
(in the limit of massless $b$-quarks) have been combined together.

The precision of theory predictions is usually quantified in terms of
renormalisation and factorisation scale variation of the NNLO results.
It is however also known that all-order effects can be sizeable and
can give rise to effects that are outside the fixed-order scale
uncertainty band. For this reason, a lot of effort is put into
combining NNLO calculations with parton shower effects, thereby
obtaining so-called NNLOPS generators. Three methods have been
suggested recently to achieve this accuracy. The {\tt UNNLOPS} approach,
which has been used for Drell-Yan~\cite{Hoeche:2014aia} and Higgs
production~\cite{Hoche:2014dla}, is based on partitioning the phase
space into an unresolved and a one-jet region and a subsequent
matching to parton shower for the resolved one-jet region.
The {\tt Geneva} approach~\cite{Alioli:2013hqa} instead uses the
next-to-next-to-leading logarithmic (NNLL) accurate resummation for a
specific observable to essentially partition the phase space. This
method has been applied recently to Drell-Yan
production~\cite{Alioli:2013hqa}.
Finally, the \noun{MiNLO}{} approach~\cite{Hamilton:2012np,Hamilton:2012rf}
relies on first using \noun{MiNLO}{} to achieve an NLO merging of the
processes with the production of the colour-singlet state ($X$) and
the same processes with one additional jet ($X+1$ jet), and on
performing a reweighing of the \noun{MiNLO}{} $X+1$ jet events to NNLO Born
distributions for $X$. This method has been applied recently to Higgs
production~\cite{Hamilton:2013fea,Hamilton:2015nsa}, Drell-Yan~\cite{Karlberg:2014qua} and \textrm{HW}{} production~\cite{Astill:2016hpa}.

In this paper we consider the production of a Higgs boson in
association with a Z boson and consider the decay of the Higgs to
bottom quarks, the decay mode with the largest branching ratio, and
the decay of the Z boson to leptons. We build a Monte Carlo that is
NNLO accurate in production, preserves NLO accuracy in the decay and
includes parton shower effects.
We also include the NNLO $gg\!\to\!\rm{HZ}{}$ channel at leading order in
production, including LO corrections in decay and parton shower
effects. This subprocesses is added separately, and we assess its
numerical impact.
Our implementation uses and adapts the \noun{POWHEG-BOX-RES}{} code, which is
based on \noun{POWHEG-BOX-V2}{} but has a resonant-aware treatment of
internal resonances~\cite{Jezo:2015aia}, and hence it is suited to
treat NLO corrections to production and decay.

The paper is organised as follows. In Sec.~\ref{sec:method} we outline
the method used to obtain \noun{NNLOPS}{} accurate predictions including the
NLO treatment of the decay. In particular we explain how the latter is
included together with \noun{MiNLO}{} within the \noun{POWHEG-BOX-RES}{}
framework. We detail how we parametrise the phase-space, and also
explain how we treat the ${\cal O}(\alpha_{\scriptscriptstyle \mathrm{s}}^2)$ $gg\!\to\!\rm{HZ}{}$ contribution. In
Sec.~\ref{sec:practical} we give details about our practical
implementation, as well as about our interface to the parton
shower. In Sec.~\ref{sec:validation} we validate our results by
checking that we reproduce NNLO results for Born-like observables,
also for variables not used for the reweighting. In
Sec.~\ref{sec:pheno} we present distributions with fiducial cuts
inspired by the recent ATLAS analysis of
ref.~\cite{Aaboud:2017xsd}. Finally, we present our conclusions in
Sec.~\ref{sec:conclu}. A number of Appendices provide more details
about the treatment of the decay at NLO, the spectral decomposition
that we use to parametrise the phase space, the dependence of the
matrix element on the extra polar angle used in the phase-space
parametrisation, and the impact of $gg\!\to\!\rm{HZ}{}$ contribution.

\section{Outline of the method}
\label{sec:method}

In this work we consider the production of a Higgs boson in
association with a Z boson, followed by the Z boson decay into a pair
of leptons and the Higgs boson decay into pair of $b$-quarks
\begin{equation}
  \label{eq:process}
  pp \longrightarrow HZ
  \longrightarrow
  \left(b\bar{b}\right)
  \left(\ell^{+}\ell^{-}\right).
\end{equation}
The decay of the Z boson is treated exactly with all spin correlations
between initial and final-state fermions taken into account. The decay
of the Higgs boson is treated in the narrow-width approximation at
next-to-leading order in QCD.

In order to achieve \noun{NNLOPS}{} accuracy we follow the method of
reweighting Les Houches events (\noun{LHE}{}), produced by \noun{MiNLO}{} improved
\scalebox{1.0}{\texttt{HZJ}}{} generator (\noun{HZJ-MiNLO}{}), with NNLO accurate fixed-order
predictions, differential in the Born phase space.  The procedure was
first proposed in~\cite{Hamilton:2012rf} and later implemented for
various colour-singlet production
processes~\cite{Hamilton:2013fea,Karlberg:2014qua,Hamilton:2015nsa,Astill:2016hpa}.
In its simplest implementation, the method consists in reweighting
\noun{LHE}{} event samples obtained with the \noun{HZJ-MiNLO}{} generator using
multi-differential \noun{HZNNLO}{} distributions, with the factor
\begin{equation}
  \label{eq:rwgt-k-factor}
  \mathcal{W}(\Phi_{\sss \ell\bar\ell b\bar{b}}{})
  =
  \frac{ d\sigma_{\scalebox{0.75}{\texttt{NNLO}}{}}(\Phi_{\sss \ell\bar\ell b\bar{b}}) }{ d\sigma_{\scalebox{0.75}{\texttt{MiNLO}}{}}(\Phi_{\sss \ell\bar\ell b\bar{b}}) },
\end{equation}
where $\Phi_{\sss \ell\bar\ell b\bar{b}}{}$ is the Born phase-space of
process~\eqref{eq:process}. The resulting event sample (which we refer
to by \noun{HZ-NNLOPS(LHEF)}{}) is NNLO accurate: by construction the method
provides NNLO accuracy for all Born distributions and 1-jet
observables remain NLO accurate since the reweighting factor differs
from one by ${\cal O}(\alpha_s^2)$ corrections.
For the proof we refer the reader
to~\cite{Hamilton:2013fea,Karlberg:2014qua}. Furthermore, a subsequent
parton shower will not spoil the claimed accuracy provided that the
hardest real radiation for each event is generated by \noun{POWHEG}{}
itself.
This procedure was applied recently to \textrm{HW}{} production in
ref.~\cite{Astill:2016hpa}, hence we refer the reader to that paper
for further details.
Instead, in this section, we first give a detailed description of the
treatment of the NLO ${\textrm{H}}\!\to\!b\bar{b}{}$ decay, which is the new element of this
work (Sec.~\ref{sec:hbb-nlo-dec}). We then give some technical details
on how the reweighting to NNLO was achieved
(Sec.~\ref{sec:phase-space-decomp}
and~\ref{sec:reweighting-procedure}), and, finally, in
Sec.~\ref{sec:gghz-treatment} we discuss the inclusion of the
loop-induced $gg\!\to\!\rm{HZ}{}$ process, which is part of the NNLO corrections
to $pp\to\textrm{HZ}{}$.

\subsection{\noun{HZJ-MiNLO}{} with ${\textrm{H}}\!\to\!b\bar{b}{}$  decay at NLO}
\label{sec:hbb-nlo-dec}

In this work we use the \noun{MiNLO}{} prescription only for the production
part of the process, and match this to a ``resonance improved''
\noun{POWHEG}{} implementation of the NLO QCD calculation of the ${\textrm{H}}\!\to\!b\bar{b}{}$
decay.
As described in the previous subsection, we treat the Higgs boson
decay in the narrow-width approximation (NWA).

We start by introducing the $\bar B$ function~\cite{Frixione:2007vw} that
we use in our \scalebox{1.0}{\texttt{HZJ}}{} code, which reads schematically
\begin{equation}
  \begin{split}
    \bar B &=
    \alpha_s(q_t^2) \Delta_q^2(Q,q_t)
    \frac{{
         {\rm{Br}}({\rm{H}}\!\rightarrow\! b\bar{b}) {} }}{\Gamma_{\scalebox{0.5}{\text{NLO}}}}
    \left[
      B_{\scriptscriptstyle  \rm HZJ} (1 - 2 \Delta_q^{(1)}(Q,q_t))
      \dgam{0} \right.
      \\ &
      + \left(V_{\scriptscriptstyle  \rm HZJ} +\int d
      \phi_r R_{\scriptscriptstyle  \rm HZJ}(\phi_r) \right) \dgam{0}
      + \left.  B_{\scriptscriptstyle  \rm HZJ} \left(\dgam{V}
      + \int d\phi_r \dgam{R}(\phi_r)\right) \right]\,,
\end{split}
\label{eq:minlo1}
\end{equation}
where the Higgs propagator is left implicit; $ {\rm{Br}}({\rm{H}}\!\rightarrow\! b\bar{b}) {}$ is the best
prediction for the Standard Model ${\textrm{H}}\!\to\!b\bar{b}{}$ branching ratio; and
$\dgam{0/V/R}$ are the Born, virtual, and real squared amplitudes for
the ${\textrm{H}}\!\to\!b\bar{b}{}$ decay, differential in their kinematics.
$\Gamma_{\scalebox{0.5}{\text{NLO}}{}} = \gam{0} + \gam{1}$ denotes the NLO accurate
${\textrm{H}}\!\to\!b\bar{b}{}$ partial decay width.
With $\Delta_q(Q,q_t)$ we denote the \noun{MiNLO}{} Sudakov form factor for
quark induced boson production (see ref.~\cite{Hamilton:2012rf} for
its precise definition) and $\Delta_q^{(1)}(Q,q_t)$ is its ${\cal
  O}(\alpha_s)$ expansion.
The hard scale in the Sudakov is set to $Q^2 = (p_Z+p_H)^2$ and $q_t$
is the transverse momentum of the $\textrm{HZ}{}$ system, where the Higgs is
obtained from the sum of the momenta of its decay products ($b \bar b$
or $b \bar b g$).
In the formula above the additional $\alpha_s$ factor in the NLO
correction is contained implicitly in the $V$ and $R$ functions, as
well as in $\dgam{V}$ and $\dgam{R}$.
In the former two, following the \noun{MiNLO}{} prescription, we set the
central renormalisation scale to $\mu_R = q_t$, whereas for the decay
we set the central scale to $\mu_{\scriptscriptstyle \mathrm{R}}= M_{\rm{H}}$ since this is the natural
scale for the decay and no \noun{MiNLO}{} procedure is applied to it (in
App.~\ref{App:higgs-nlo-decay} we denote as $\mu_r$ the
renormalisation scale for the decay).

If we integrate Eq.~\eqref{eq:minlo1} over the phase space of all final-state light
partons we obtain
\begin{eqnarray}
 \label{eq:MiNLO-hz}
 \hspace*{-1.0cm}d\sigma_{\scalebox{0.75}{\texttt{MiNLO}}}(\Phi_{\sss \ell\bar\ell b\bar{b}}{}) &=&  {\rm{Br}}({\rm{H}}\!\rightarrow\! b\bar{b})  \!\cdot
 \!\left[ \left( \dsig{0}_{\rm \scriptscriptstyle  HZ} + \dsig{1}_{\rm \scriptscriptstyle
     HZ}\right) \cdot\frac{ \dgam{0} +
     \dgam{1}}{\Gamma_{\scalebox{0.5}{\text{NLO}}}} \!+
   d\tilde{\sigma}^{(2)}_{\rm \scriptscriptstyle  HZ} \cdot \frac{\dgam{0}}{\Gamma_{\scalebox{0.5}{\text{NLO}}}} \right]+{\cal O}(\alpha_s^3),
\end{eqnarray}
where
\begin{equation}
  \qquad
  \dgam{1}
  =
  \dgam{V}
  + \int d\phi_r \dgam{R}(\phi_r)\,,\,
\end{equation}
and $d \sigma_{\rm \scriptscriptstyle  HZ}^{(i)}$ denotes the ${\cal O}(\alpha_s^i)$
correction to the \textrm{HZ}{} production cross section.  The
$d\tilde\sigma^{(2)}$ denotes the ${\cal O}(\alpha_s^2)$ part of the
\noun{HZJ-MiNLO}{} computation, which corresponds to double-real and
real-virtual parts of $\textrm{HZ}{}$ production at NNLO.

We obtain NNLO prediction (\emph{without} the loop-induced $gg\!\to\!\rm{HZ}{}$
contribution, which is discussed in Sec.~\ref{sec:gghz-treatment}) for
the production combined with NLO corrections to the decay from
\texttt{MCFM-8.0}{}, whose output is
\begin{eqnarray}
  \label{eq:NNLO-mcfm}
  d\sigma_{\scalebox{0.75}{\texttt{NNLO}}{}}(\Phi_{\sss \ell\bar\ell b\bar{b}}{}) &=&  {\rm{Br}}({\rm{H}}\!\rightarrow\! b\bar{b}) {} \cdot \left[
    \dsig{0}_{\rm \scriptscriptstyle  HZ} \cdot \frac{\dgam{0} +
      \dgam{1}}{\Gamma_{\scalebox{0.5}{\text{NLO}}}} + (\dsig{1}_{\rm \scriptscriptstyle
      HZ} + \,\dsig{2}_{\rm \scriptscriptstyle  HZ}) \cdot \frac{\dgam{0}}{\gam{0}} \right]\,.
\end{eqnarray}

It is easy to check that after integrating out the decay of the Higgs
boson in equation~\eqref{eq:NNLO-mcfm} one
recovers the fully inclusive NNLO result multiplied by the overall branching
ratio. One can also easily verify that
\begin{equation}
  \mathcal{W}(\Phi_{\sss \ell\bar\ell b\bar{b}}{}) = \frac{
    d\sigma_{\scalebox{0.75}{\texttt{NNLO}}{}}(\Phi_{\sss \ell\bar\ell b\bar{b}}{}) }{
    d\sigma_{\scalebox{0.75}{\texttt{MiNLO}}}(\Phi_{\sss \ell\bar\ell b\bar{b}}{}) } = 1 + \frac{\left(\sigma^{(2)}
    - \tilde{\sigma}^{(2)}\right) }{\sigma^{(0)}} +
  \mathcal{O}\left(\alpha_{\scriptscriptstyle \mathrm{s}}^{3}\right)\,,
\end{equation}
which means that reweighting \noun{HZJ-MiNLO}{} events with this factor does
not spoil the NLO accuracy of the event sample in the $\textrm{HZ+jet}{}$
region, since the rescaling is equal to one up to $\mathcal{O}(\alpha_{\scriptscriptstyle \mathrm{s}}^2)$
terms. In the following section we describe how we proceed to obtain
distributions differential in the Born phase space $\Phi_{\sss \ell\bar\ell b\bar{b}}{}$.

\subsection{Phase-space parametrisation}
\label{sec:phase-space-decomp}
Our Born phase space contains four final-state particles, two leptons
($\ell^{+}, \ell^{-}$) and two bottom quarks ($b, \bar b$).
After neglecting an irrelevant azimuthal angle and upon inclusion of
the initial-state degrees of freedom we are left with 9 independent
dimensions. Furthermore we can factorise the Born phase-space as
follows:
\begin{equation}
   d\Phi_{\sss \ell\bar\ell b\bar{b}} = d\Phi_{\sss\rm{H}\ell\bar\ell} \times (2\pi)^3 dq^2\times d\Phi_{\sss H \to b\bar{b}},
\end{equation}
where $q^2$ is the virtuality of the Higgs boson, $\Phi_{\sss\rm{H}\ell\bar\ell}$ is the
6-dimensional phase space for the production of an undecayed Higgs
boson with a pair of leptons from the decay of the associated Z boson,
and $\Phi_{\sss H \to b\bar{b}}$ is the 2-dimensional phase space for the decay of a
Higgs boson into a pair of $b$-quarks.  By working in the NWA we
perform a reweighting only on the first part of the phase-space
$\Phi_{\sss\rm{H}\ell\bar\ell}$, as the decay is already treated at the required accuracy
(NLO) in our implementation of \noun{HZJ-MiNLO}{}.

A parametrisation of $\Phi_{\sss\rm{H}\ell\bar\ell}{}$ can be defined in a number of
ways. After careful considerations we have chosen the invariant mass
and rapidity of the \textrm{HZ}{} system ($M_{\sss \HZ}$ and $y_{\sss \HZ}$) to be the first
two variables. As a third variable we choose $\cos\!\alpha$, where $\alpha$
is the polar angle of the Z boson with respect to the beam axis in the
frame where the \textrm{HZ}{} system is at rest,~\emph{i.e.}
\begin{equation}
  \label{eq:cosx-def}
  \cos\!\alpha = \frac{ \vec{p}\,'_{Z} \cdot \hat{z}' }{ |\vec{p}\,'_Z|\; |\hat{z}'| },
\end{equation}
where $'$ indicates that directions are expressed in the \textrm{HZ}{}
rest-frame and the original $\hat{z}$ direction is along the beam
axis. Subsequently we choose the invariant mass of the Z boson,
$M_{\sss \ell\bar\ell}{}$, and a convenient choice for the last two dimensions is to use
Collins-Soper angles ($\theta^*,\phi^*$) defined
in~\cite{Collins:1977iv}. In summary the full phase space parametrisation reads
\begin{equation}
  \label{eq:ps-coord}
  \Phi_{\sss\rm{H}\ell\bar\ell} = \left\{ M_{\sss \HZ}, y_{\sss \HZ}, \cos\!\alpha, M_{\sss \ell\bar\ell}{}, \theta^*, \phi^* \right\}.
\end{equation}

Following the arguments of~\cite{Collins:1977iv} and the discussion in
Sec.~2 of~\cite{Astill:2016hpa} we
parametrise the ($\theta^*,\phi^*$)-dependence as follows:
\begin{eqnarray}
  \label{eq:diff-xs}
  \frac{d\sigma}{d\Phi_{\sss\rm{H}\ell\bar\ell}} &=& \frac{d^6\sigma}{dM_{\sss \HZ}{}\, dy_{\sss \HZ}{}\,
    d\cos\!\alpha\, d\cos \theta^* d\phi^*} \nonumber \\ &=& \frac{3}{16\pi}
  \left ( \frac{d \sigma}{d\Phi_{\sss \HZ}}(1+\cos^2\theta^*) +
  \sum_{i=0}^{7} A_i(\Phi_{\sss \HZ} ) f_i(\theta^*, \phi^*) \right)\,,
\end{eqnarray}
where in the second line we have used a short notation for the phase-space
without leptonic angles
\begin{equation}
  \label{eq:phase-space-3}
  \Phi_{\sss \HZ} = \left\{ M_{\sss \HZ}, y_{\sss \HZ}, \cos\!\alpha \right\}.
\end{equation}
The complete set of functions $f_{i}(\theta^*,\phi^*)$, together with
a procedure for extracting coefficients $A_{i}(\Phi_{\sss \HZ})$, is given
in equations (2.3-2.4) of ref.~\cite{Astill:2016hpa}. We note that for
practical purposes we will use only the first five $A_i$ coefficients
($A_0, \ldots A_4$) and, for simplicity, we neglect the remaining
three ($A_5,A_6,A_7$) since their impact on any distribution that we
examined is numerically negligible.

Finally, we can parametrise the dependence on the Z boson polar angle
$\alpha$ (see Eq.~\eqref{eq:cosx-def}) using a set of orthonormal
functions $g_{j}(\cos\!\alpha)$. The definition of the basis elements $g_j$
is given in App.~\ref{App:GramSchmidt}. With this choice we have
\begin{eqnarray}
  \frac{d \sigma}{d\Phi_{\sss \HZ}} = \sum_{j=0}^{N} c_j(\Phi) \,
  g_j\left(\cos\!\alpha\right), \notag\\ A_{i}\left(\Phi_{\sss \HZ}\right) =
  \sum_{j=0}^{N} a_{ij}(\Phi) \, g_j\left(\cos\!\alpha\right),
  \label{eq:poly-decomp}
\end{eqnarray}
where $c_j$ and $a_{ij}$ are coefficients depending on $\Phi=\left\{
M_{\sss \HZ}, y_{\sss \HZ} \right\}$ and $N$ is the upper limit of the sum, which in
general can be inferred by analysing the matrix elements contributing
to the cross section.  We investigate the matrix elements in
App.~\ref{App:HadronicTensor} and find that, at most, polynomials of
5th-degree in the $\cos\!\alpha$ and $\sin{\!\alpha}$ variables can
appear. Accordingly, we set $N=10$ in Eq.~\eqref{eq:poly-decomp}.

To summarise, our parametrisation of the fully differential cross
section as used in the reweighting procedure reads
\begin{equation}
  \label{eq:diff-reconstruction}
  \frac{d\sigma}{d\Phi_{\sss\rm{H}\ell\bar\ell}} = \frac{3}{16\pi} \sum_{j=0}^{10}\left (
  c_j(\Phi)\,(1+\cos^2\theta^*) + \sum_{i=0}^{7} a_{ij}(\Phi)\,
  f_i(\theta^*, \phi^*) \right)g_{j}\left(\cos\!\alpha\right)\,,
\end{equation}
where the functions $g_j(\cos\!\alpha)$ are defined in Eq.~\eqref{eq:gdef}
and the coefficients $c_j(\Phi)$ and $a_{ij}(\Phi)$ can be obtained
from Eq.~\eqref{eq:poly-decomp} by exploiting the orthonormality of
the $g_j$ functions.

\subsection{Reweighting procedure}
\label{sec:reweighting-procedure}
The reweighting procedure leaves some degree of freedom for the
user. The simple rescaling with a factor, presented
in~\eqref{eq:rwgt-k-factor}, spreads the corrections uniformly over
the whole phase-space. However we know that regions where the
\textrm{HZ}{} system is accompanied by hard QCD radiation is equally well
described by both predictions, \noun{HZNNLO}{} and \noun{HZJ-MiNLO}{}. Hence, it is
desirable to limit the corrections to the phase space with no hard
jet.
To achieve this, we proceed along the lines of the prescription
presented in~\cite{Hamilton:2013fea}. We split the cross section into
two parts
\begin{equation}
  \label{eq:sigmaA-sigmaB}
  d\sigma_A = d\sigma\, h(p_t), \qquad d\sigma_B = d\sigma \left(1 - h(p_t) \right),
\end{equation}
where
\begin{equation}
  \label{eq:h-function}
  h(p_t) = \frac{ (M_{\rm{H}}+M_{\rm{Z}})^2 }{ (M_{\rm{H}}+M_{\rm{Z}})^2 + p_t^{\,2} }\,,
\end{equation}
and $p_t$ is the transverse momentum of the leading jet, computed here
using the $k_t$-algorithm with $R=0.4$. With such a choice
Eq.~\eqref{eq:rwgt-k-factor} takes form
\begin{eqnarray}
  \label{eq:W-rwgt}
  \mathcal{W}(\Phi_{\sss\rm{H}\ell\bar\ell}{},p_t) &=& h(p_t)\,\frac{ \int
    d\sigma_{\scalebox{0.75}{\texttt{NNLO}}{}} \,\delta(\Phi_{\sss\rm{H}\ell\bar\ell} - \Phi_{\sss\rm{H}\ell\bar\ell}(\Phi)) - \int
    d\sigma_{\scalebox{0.75}{\texttt{MiNLO}}{},B}\, \delta(\Phi_{\sss\rm{H}\ell\bar\ell} - \Phi_{\sss\rm{H}\ell\bar\ell}(\Phi)) }{
    \int d\sigma_{\scalebox{0.75}{\texttt{MiNLO}}{},A}\, \delta(\Phi_{\sss\rm{H}\ell\bar\ell} - \Phi_{\sss\rm{H}\ell\bar\ell}(\Phi)) }
  \notag\\ & & + \left(1-h(p_t)\right).
\end{eqnarray}
This procedure smoothly turns off the corrections when moving to
regions of phase space with hard emissions whilst preserving the total
cross section,
\begin{equation}
  \left(\frac{d\sigma}{d\Phi_{\sss\rm{H}\ell\bar\ell}}\right)_{\scalebox{0.75}{\texttt{NNLOPS}}{}}
  =
  \left(\frac{d\sigma}{d\Phi_{\sss\rm{H}\ell\bar\ell}}\right)_{\scalebox{0.75}{\texttt{NNLO}}{}}.
\end{equation}
It is worth noting that choosing the transverse momentum $p_t$ in
Eq.~\eqref{eq:h-function} as the transverse momentum of hardest jet is
not equivalent to choosing the transverse momentum of the
colour-singlet, the difference being due to configurations where the
colour-singlet has small transverse momentum and it is accompanied by
hard QCD emissions whose transverse momenta counterbalance each other.
The latter configurations dominate in the region where $p_{t,\textrm{HZ}{}} \sim
0$.

In the next section we will explain how we have included in our
simulation the loop-mediated $gg\!\to\!\rm{HZ}{}$ contribution, which was not
included in ${d\sigma}_{\scalebox{0.75}{\texttt{NNLO}}{}}$ (and hence not even in the
reweighting factor defined in Eq.~\eqref{eq:W-rwgt}), whilst being
formally $\mathcal{O}(\alpha_{\scriptscriptstyle \mathrm{s}}^2)$.

\subsection{Treatment of the $gg\!\to\!\rm{HZ}{}$ contribution}
\label{sec:gghz-treatment}
The $\mathcal{O}(\alpha_{\scriptscriptstyle \mathrm{s}}^2)$ contributions of the form $gg\!\to\!\rm{HZ}{}$, that
appear in the process of Higgs boson production in association with
a Z boson, originate from the squared 1-loop diagrams
shown for example in Fig.~5 of ref.~\cite{Campbell:2016jau}.
There are two classes of contributions: box diagrams involving a Yukawa
coupling of the Higgs boson, and triangle diagrams where the Higgs
boson is radiated from the Z boson which couples to the fermion
loop.
Both contributions vanish trivially when massless quarks of the first
and second generation run in the loop.
These $gg\!\to\!\rm{HZ}{}$ contributions are known to constitute a significant
part of the cross
section~\cite{Hespel:2015zea,Goncalves:2015mfa,Ferrera:2014lca,Campbell:2016jau,Ferrera:2017zex},
especially when the invariant mass of the produced \textrm{HZ}{} system is
larger than twice the top-quark mass. Being loop-induced, so far only
approximate NLO corrections are known for this
channel~\cite{Altenkamp:2012sx,Harlander:2014wda}. On the other hand,
their loop origin makes these terms particularly sensitive to New
Physics~\cite{Englert:2013vua,Harlander:2013mla,Harlander:2018yns}.

In our reweighting procedure we do not include this contribution, but
we will treat it independently using a separate event sample produced
using a leading order implementation of loop-induced $gg\!\to\!\rm{HZ}{}$ process
implemented in \noun{POWHEG}{}~\cite{Luisoni:2013kna}, which includes top
and bottom quarks in the loop.\footnote{The code can be obtained from
  \texttt{svn://powhegbox.mib.infn.it/trunk/User-Processes-V2/ggHZ}.}
We note that the $gg\!\to\!\rm{HZ}{}$ contribution can be treated separately
since it is finite.  Furthermore, parton shower radiation from
gluon-initiated hard processes is typically different from processes
also involving initial-state quarks. From that point of view, it is
important not to include the $gg\!\to\!\rm{HZ}{}$ contribution through a simple
reweighting.
Further discussion of the $gg\!\to\!\rm{HZ}{}$ channel is presented in
App.~\ref{App:gghz}.
For this contribution we do not include any radiative correction to
the the ${\textrm{H}}\!\to\!b\bar{b}{}$ decay, hence the radiation from the decay is
fully taken care of by \PYTHIA{8}. Higher-order NLO corrections to this
decay could also be included with relatively little effort.

\section{Practical implementation}
\label{sec:practical}

In this section we first discuss the codes used to obtain our
predictions, as well as the relevant settings and the parameters. We
also describe subtleties related to the interface to \PYTHIA{8} when
radiation from resonances is taken into account.

\subsection{Codes and settings}
\label{sec:setting}

In order to obtain an ensemble of NNLOPS accurate Les Houches events
for the process in Eq.~\eqref{eq:process} we need fully differential
predictions from an NNLO fixed-order calculation, and an NLO accurate
event-sample for \scalebox{1.0}{\texttt{HZJ}}{} production improved with the \noun{MiNLO}{}
prescription.

For the NNLO fixed-order prediction we use the \texttt{MCFM-8.0}{}
calculation~\cite{Campbell:2016jau}. In order to obtain both the NNLO
accuracy for the production of the \textrm{HZ}{} resonance as well as the NLO
accuracy of hadronic decay of Higgs boson, ${\textrm{H}}\!\to\!b\bar{b}{}$, we
perform two separate runs of the program with \texttt{nproc=101} at
the \texttt{'nnlo'} order (for the first) and \texttt{nproc=1010} at
the \texttt{'nlo'} order (for the latter). The prediction presented in
Eq.~\eqref{eq:NNLO-mcfm} is simply obtained by adding results of the two runs.
As it was pointed out in Sec.~\ref{sec:hbb-nlo-dec} and
Sec.~\ref{sec:gghz-treatment}, we do not include $gg\!\to\!\rm{HZ}{}$
contributions in the reweighting procedure. To remove them we have
modified part of the \texttt{MCFM-8.0}{} code, which computes double-virtual
corrections. We will include this contribution in our phenomenological
analysis in Sec.~\ref{sec:pheno}, as stated clearly in the appropriate
places.

The initial sample of Les Houches events is generated using the
\noun{HZJ-MiNLO}{} package, originally in \noun{POWHEG-BOX-V2}{}, which we adapted
to run in \noun{POWHEG-BOX-RES}{}~\cite{Jezo:2015aia}.
We have also extended the original package to include NLO corrections
of the Higgs boson decay into pair of $b$-quarks, as discussed in
Sec.~\ref{sec:hbb-nlo-dec}. The relevant matrix elements have been
reported in App.~\ref{App:higgs-nlo-decay}.
Despite the fact that there is no interference between production and
decay, in order to treat the radiation from the resonance we have made
use of the \noun{POWHEG-BOX-RES}{} framework~\cite{Jezo:2015aia}.

In our work we use \texttt{PDF4LHC15\_nnlo\_mc} parton distribution
functions~\cite{Ball:2014uwa,Harland-Lang:2014zoa,Dulat:2015mca,Carrazza:2015hva}. We
set $M_{\rm{H}} = 125.0$ GeV, $\Gamma_H = 4.088$ MeV, $M_{\rm{Z}}=91.1876$ GeV, and
$\Gamma_Z = 2.4952$ GeV. Moreover $G_F = 1.16639\cdot 10^{-5}
\rm{GeV}^{-2}$, $\sin^2\theta_{W} = 0.2223$,
$\alpha_{\scalebox{0.5}{\rm{EM}}}(M_{\rm{Z}}) = 128.89$, and $ {\rm{Br}}({\rm{H}}\!\rightarrow\! b\bar{b}) {} =
0.5824$.

For the contributions where the Higgs boson is radiated from a
heavy-quark loop we use pole mass of the heavy quark as it is usually
done in publicly available codes~\cite{Campbell:2016jau}. In
particular we set pole mass of the bottom quark to $m_b = 4.92$ GeV
and pole mass of the top quark to $m_t = 173.2$ GeV. Moreover, for the
bottom Yukawa coupling in ${\textrm{H}}\!\to\!b\bar{b}{}$ decay we use its $\overline{\rm MS}$ running
mass evaluated at scale $M_{\rm{H}}$. The running masses are computed from
the pole masses using an ${\cal O}(\alpha_s^2)$
conversion~\cite{Chetyrkin:2000yt} and the numerical value of the
$\overline{\rm MS}$ mass for the bottom quark, obtained by running the strong
coupling using \noun{LHAPDF}{}, is $m_b(M_{\rm{H}}) = 3.16$ GeV.

The NNLO fixed-order prediction is obtained using fixed
renormalisation and factorisation scale equal to sum of the Higgs
boson and the Z boson mass, $\mu = M_{\rm{H}}+M_{\rm{Z}}$. The scale choice in
\noun{HZJ-MiNLO}{} case for the production is fixed by \noun{MiNLO}{}
procedure~\cite{Hamilton:2012np}, while for the decay the central
scale choice is $M_{\rm{H}}$, as explained in Sec.~\ref{sec:hbb-nlo-dec}.
We estimate the theoretical uncertainty using 7 point scale variation
for both fixed-order NNLO results as well as \noun{MiNLO}{} predictions. The
scale variation in production and decay are always correlated (this
includes the decay renormalisation scale,~\emph{i.e.} the scale at which we
evaluate the $\overline{\rm MS}$ $b$-quark mass).
We perform correlated variations in \texttt{MCFM-8.0}{} and
\noun{HZJ-MiNLO}{}, that is, our final uncertainty is an envelop of 7
scale combinations,~\emph{i.e.} for a given $(K_r,K_f)$ choice in \noun{HZJ-MiNLO}{}
we use the same choice in \texttt{MCFM-8.0}{} and as usual we consider variations
of the central scale by a factor two up and down, restricted to $1/2
\le K_r/K_f < 2$.
The $gg\!\to\!\rm{HZ}{}$ contribution is then added with the same $(K_r,K_f)$
choice.

When interfacing our fixed-order predictions to a parton shower we use
\PYTHIA{8}~\cite{Sjostrand:2014zea}, as detailed more precisely in the
next subsection. Unless stated otherwise, predictions are shown at
parton level, with no multi-parton interactions.

\subsection{Interface to Parton Shower}
\label{sec:intshower}
In order to combine our results with a parton shower we follow an
approach similar to the one first introduced for the NLO $t\bar t$
production treatment in \cite{Campbell:2014kua}, that allows for a
generation of radiation also from resonances. In our simulation we set
the flag \texttt{allrad} to 0, which means that the NLO \noun{POWHEG}{}
emission is generated at most from one singular region, associated
either with the production stage or with the radiation from a
resonance and we do not consider radiation from multiple regions.
This is the standard \noun{POWHEG}{} procedure to generate the hardest
radiation. In this configuration \noun{POWHEG}{} uses the usual highest bid
mechanism to choose the origin of the emission.

For the parton shower we use \PYTHIA{8}.  A requirement for the
matching to the parton shower to work properly is that the hardest
radiation should be the one generated by \noun{POWHEG}{}. This is usually
achieved by setting a value of \texttt{scalup} in every event, which
sets the veto scale for the parton shower. One subtlety is however
that the definition of the hardness of the radiation from the decay in
\noun{POWHEG}{} and \PYTHIA{8} differ. As a consequence, after an event is
showered, we recompute the hardness of the first emission generated by
\PYTHIA{8} using the \noun{POWHEG}{} formula and accept the showered event
only if this hardness is lower than the \texttt{scalup} value of the
given event. If this is not the case, we shower the event again until
the new showered event meets the required condition. Details of the
hardness definition used in \noun{POWHEG}{} and \PYTHIA{8} are given in
App.~A of ref.~\cite{Campbell:2014kua}. We have checked that our
procedure to veto radiation gives results that are fully compatible
with those obtained by the procedure encoded in the
\noun{PowhegHooks}{}-class provided by \PYTHIA{8}~\cite{Sjostrand:2014zea}.

\section{Validation}
\label{sec:validation}

In the following section we present the validation of our results.  We
carefully compare distributions prepared from reweighted Les Houches
events (\noun{HZ-NNLOPS(LHEF)}{}) with the ones obtained using fixed-order NNLO
code (\texttt{MCFM-8.0}{}). We remind the reader that, in the reweighting
procedure, we don't take into account one-loop squared contributions
arising from $gg\!\to\!\rm{HZ}{}$ channel, as specified and motivated in
Sec.~\ref{sec:reweighting-procedure} and~\ref{sec:gghz-treatment}
respectively.  Therefore all the plots of this section do not contain
the $gg\!\to\!\rm{HZ}{}$ channel, which instead will be included in
Sec.~\ref{sec:pheno}.

In the plots of this section, the blue, green and red markers
represent results from \noun{HZJ-MiNLO}{} Les Houches events, the fixed-order
calculation obtained with \texttt{MCFM-8.0}{}, and the reweighted \noun{HZ-NNLOPS(LHEF)}{} event
sample, respectively. The uncertainty band represents the usual scale
variation uncertainty, as described in detail in the previous section.

The first pair of plots that we want to present is the distribution of
the invariant mass of final-state leptons $M_{\sss \ell\bar\ell}$ and the distribution
of the rapidity of the \textrm{HZ}{} system $y_{\sss \HZ}$, which are shown in the
left and right panel of Fig.~\ref{fig:mll-yHZ-plots}, respectively.
\begin{figure}
  \centering
  \includegraphics[page=1,width=205pt]{graphics/validation/born-mll.pdf}
  \hspace*{0.25cm}
  \includegraphics[page=1,width=205pt]{graphics/validation/born-yhz.pdf}
  \caption{The differential distributions of the invariant mass of
    final-state leptons $M_{\sss \ell\bar\ell}$ (left panel) and the distribution of
    the rapidity of the \textrm{HZ}{} system $y_{\sss \HZ}$ (right panel). The one-loop squared terms
    from the $gg\!\to\!\rm{HZ}{}$ channel have not been included.}
  \label{fig:mll-yHZ-plots}
\end{figure}
We start by noting that the ratio of \texttt{MCFM-8.0}{} to \noun{HZJ-MiNLO}{}, bottom
panel of Fig.~\ref{fig:mll-yHZ-plots}, is constant, which is along the
lines of our assumption that the reweighting factor
$\mathcal{W}(\Phi_{\sss \ell\bar\ell b\bar{b}}{})$ should be constant along the $M_{\sss \ell\bar\ell}$
direction. We do not repeat the thorough procedure of validation, which
was included in our previous work~\cite{Astill:2016hpa}.
The distribution of $y_{\sss \HZ}$, right panel of
Fig.~\ref{fig:mll-yHZ-plots}, is again properly reproduced by our
calculation across the whole spectrum. We take note of the fact that
the discrepancies at the edges of the distribution are in the regions
of phase-space which are poorly populated. More precisely, having used
distributions with varying bin-size for the reweighting, all events
with $y_{\sss \HZ} \lesssim -3$ (or $y_{\sss \HZ} \gtrsim +3$) fall into the first (or
the last) bin of the differential reweighting factor
$\mathcal{W}(\Phi_{\sss \ell\bar\ell b\bar{b}}{})$.
The description of the forward rapidity region can be improved by
increasing the statistics of the multi-differential distributions and
by including more bins at large rapidity.
In Fig.~\ref{fig:mHZ-plots}, we present the differential distribution
of $M_{\sss \HZ}$ in two different mass regions.
\begin{figure}
  \centering
  \includegraphics[page=1,width=205pt]{graphics/validation/born-mhz.pdf}
  \hspace*{0.25cm}
  \includegraphics[page=1,width=205pt]{graphics/validation/born-mhz-zoom.pdf}
  \caption{The differential distributions of the invariant mass of the
    \textrm{HZ}{} system $M_{\sss \HZ}$ in two different mass regions. The one-loop
    squared terms from $gg\!\to\!\rm{HZ}{}$ channel have not been included.}
  \label{fig:mHZ-plots}
\end{figure}
For this distribution the
difference between \noun{HZJ-MiNLO}{} and the NNLO is small and flat over the
whole range. After reweighting, we find perfect agreement between NNLO
and \noun{HZ-NNLOPS(LHEF)}{} results.

As the next step, we look closely at the differential distributions of
the angular variables: $\cos\!\alpha$ and Collins-Soper angles. The
distribution of $\cos\!\alpha$ is presented in Fig.~\ref{fig:cosx-plots}. We
recollect that this dependence was not just recorded as a histogram,
but rather parametrised in terms of spectral modes,
Eq.~\eqref{eq:poly-decomp}. This has improved the stability of the
distribution, and as a consequence of the reweighting factor, which is
very useful when working with samples of limited statistics.
\begin{figure}
  \centering
  \includegraphics[page=1,width=205pt]{graphics/validation/born-cosx.pdf}
  \hspace*{0.25cm}
  \includegraphics[page=1,width=205pt]{graphics/validation/born-a2.pdf}
  \caption{The differential distributions of the Z boson polar angle
    with respect to the beam axis, defined in Eq.~\eqref{eq:cosx-def}:
    differential cross section as a function of $\cos\!\alpha$ (left panel)
    and the dependence of coefficient $A_2(\cos\!\alpha)$, see
    Eq.~\eqref{eq:diff-xs}. The one-loop squared terms from the
    $gg\!\to\!\rm{HZ}{}$ channel have not been included.}
  \label{fig:cosx-plots}
\end{figure}
Further we check the quality of the reconstruction of the
Collins-Soper angles. We present the relevant distributions in
Fig.~\ref{fig:cs-angles}.
In summary, we find that for all variables used for the reweighting,
the NNLO and \noun{HZ-NNLOPS(LHEF)}{} predictions agree within their statistical
fluctuations.
\begin{figure}
  \centering
  \includegraphics[page=1,width=205pt]{graphics/validation/born-thCS.pdf}
  \hspace*{0.25cm}
  \includegraphics[page=1,width=205pt]{graphics/validation/born-phiCS.pdf}
  \caption{The differential distributions of Collins-Soper angles:
    $\theta^*$ (left) and $\phi^*$ (right). The one-loop squared terms
    from the $gg\!\to\!\rm{HZ}{}$ channel have not been included.}
  \label{fig:cs-angles}
\end{figure}

To complete the validation, we also need to examine Born-like
observables that were not used in the reweighting procedure. As such,
we chose to look at transverse momentum and rapidity of the Higgs
boson (see Fig.~\ref{fig:higgs-plots}). We again confirm that both the
central values and scale variation bands are properly reconstructed
within statistical fluctuations, which increase at high transverse
momentum ($p_{t,\rm{H}} \gtrsim 400$ GeV) and large rapidity ($|y_{\rm
  H}| \gtrsim 3$).
The agreement between NNLO and \noun{HZ-NNLOPS(LHEF)}{} in these corners of phase space
could be improved further by increasing the statistics of the
reweighting factor and decreasing the bin-sizes in this region.
\begin{figure}
  \centering
  \includegraphics[page=1,width=205pt]{graphics/validation/born-pth.pdf}
  \hspace*{0.25cm}
  \includegraphics[page=1,width=205pt]{graphics/validation/born-yh.pdf}
  \caption{The differential distributions of the transverse momentum
    (left panel) and the rapidity (right panel) of the Higgs
    boson. The one-loop squared terms
    from the $gg\!\to\!\rm{HZ}{}$ channel have not been included.}
  \label{fig:higgs-plots}
\end{figure}

Finally we turn to the discussion of the distribution of the
transverse momentum of \textrm{HZ}{} system, an observable which is singular
at Born level but receives corrections due to QCD radiation at
higher-orders in perturbation theory.
We compare results obtained using two different reweighting
prescriptions: the one described in
Sec.~\ref{sec:reweighting-procedure}, presented in the left plot of
Fig.~\ref{fig:pthz-plot}, and a setup where we set the function
$h(p_t)\equiv 1$ in Eqs.\eqref{eq:sigmaA-sigmaB}-\eqref{eq:W-rwgt},
shown in the right hand side of Fig.~\ref{fig:pthz-plot}.
As expected, we observe that the \noun{HZJ-MiNLO}{} and \noun{HZ-NNLOPS(LHEF)}{}
predictions feature a Sudakov damping at low transverse momentum,
while the NNLO prediction diverges in this region.
Furthermore, we observe that for the ${h(p_t)=1}$ case, the
\noun{HZ-NNLOPS(LHEF)}{} results are uniformly shifted with respect to the
original event sample \noun{HZJ-MiNLO}{}, as the reweighting factor
$\mathcal{W}(\Phi_{\sss\rm{H}\ell\bar\ell}{})$ does not take into account any QCD
radiation. Instead, when the reweighting factor depends on the
transverse momentum of the leading jet, \noun{HZ-NNLOPS(LHEF)}{} approaches the
\noun{HZJ-MiNLO}{} curve at high-$p_{\rm t,HZ}$ values, as the effects of the
reweighting are concentrated in the region of phase-space close to the
Born kinematics, the natural habitat of $\mathcal{O}(\alpha_{\scriptscriptstyle \mathrm{s}}^2)$ virtual
corrections.
In this case, the \noun{HZ-NNLOPS(LHEF)}{} prediction at high transverse momentum
agrees with the \noun{HZJ-MiNLO}{} prediction, rather than with the pure NNLO
result. We note that in this region, all predictions are only NLO
accurate and that the former has a dynamical scale, dictated by the
\noun{MiNLO}{} prescription, while the NNLO uses a fixed renormalisation and
factorisation scale choice, $M_{\rm{H}}+M_{\rm{Z}}$.
Comparing the middle panels of Fig.~\ref{fig:pthz-plot}, it might seem
that the choice of a uniform reweighting provides a better description
of the hard part of $p_{t,\textrm{HZ}{}}$ distribution, but the apparent
agreement between \texttt{MCFM-8.0}{} and \noun{HZ-NNLOPS(LHEF)}{} results around 400-500 GeV
is accidental. In fact, at even higher transverse momenta the NNLO
result is above our \noun{HZ-NNLOPS(LHEF)}{} prediction. This behaviour is
entirely due to the aforementioned difference in scale choice.
\begin{figure}
  \centering
  \label{fig:pthz-plot}
  \includegraphics[page=1,width=205pt]{graphics/validation/born-pthz.pdf}
  \hspace*{0.25cm}
  \includegraphics[page=1,width=205pt]{graphics/validation/born-pthz-h1.pdf}
  \caption{The differential distribution of the transverse momentum of
    the \textrm{HZ}{} system results when reweighting with damping factor
    $h(p_t)$ (left panel) or without (right panel). The one-loop squared terms
    from the $gg\!\to\!\rm{HZ}{}$ channel have not been included.}
  \label{fig:pthz-plot}
\end{figure}

\section{Phenomenological results}
\label{sec:pheno}

In this section we turn to the discussion of the phenomenological
results obtained with our new code. We stress again that we consider
the production of a Higgs boson in association with a Z boson and
their subsequent decays ${\textrm{H}}\!\to\!b\bar{b}{}$ and $\rm{Z}\to
\ell^{+}\ell^{-}$, where in the following $\ell$ denotes a single
leptonic species,~\emph{e.g.} $e$ or $\mu$.  For the Higgs boson decay we
include NLO QCD corrections.

We consider 13 TeV LHC collisions.  We consider fiducial cuts inspired
by the recent ATLAS analysis of ref.~\cite{Aaboud:2017xsd}. We require
two charged leptons with $|y_{\ell}| < 2.5$ and $p_{t,\ell} > $ 7 GeV,
moreover the harder lepton should satisfy $p_{t,\ell} > $ 27 GeV. We
impose that the invariant mass of the leptons satisfies the condition
81 GeV $ < M_{\sss \ell\bar\ell} < $ 101 GeV.  Additionally we require at least
two $b$-jets with $|\eta_{j}| < 2.5$ and $p_{t,j} > $ 20 GeV.  Unless
stated otherwise, jets are defined using the flavour-$k_t$
algorithm~\cite{Banfi:2006hf} with $R=0.4$. In the flavour-$k_t$
algorithm we only consider $b$-quarks to be flavoured, and all other
light quarks to be flavourless. Using $b$-tagging, such an algorithm
can be implemented in experimental analyses.
\begin{table}
\centering
\begin{tabular}{|l|c|c|c|c|}
  \hline
\bf{Fiducial cross section}  &  \noun{HZJ-MiNLO}{} & \texttt{MCFM-8.0}{} & \noun{HZ-NNLOPS(LHEF)}{} & \noun{HZNNLOPS}{} \\
\hline
no $gg\!\to\!\rm{HZ}{}$
& $6.59^{+7.2\%}_{-6.2\%}$ fb
& $7.14^{+0.5\%}_{-0.9\%}$ fb
& $7.14^{+0.3\%}_{-0.4\%}$ fb
& $6.49^{+0.8\%}_{-0.6\%}$ fb
\\
\hline
with $gg\!\to\!\rm{HZ}{}$
& --
& $7.92^{+2.0\%}_{-1.5\%}$ fb
& $7.90^{+2.8\%}_{-2.0\%}$ fb
& $7.16^{+3.1\%}_{-2.1\%}$ fb
\\
\hline
no $gg\!\to\!\rm{HZ}{}$, high-$p_{t,Z}$
& $1.13^{+5.9\%}_{-5.3\%}$ fb
& $1.21^{+0.1\%}_{-0.2\%}$ fb
& $1.21^{+0.2\%}_{-0.3\%}$ fb
& $1.13^{+1.5\%}_{-1.2\%}$ fb
\\
\hline
with $gg\!\to\!\rm{HZ}{}$, high-$p_{t,Z}$
& --
& $1.49^{+5.3\%}_{-4.1\%}$ fb
& $1.48^{+5.3\%}_{-4.0\%}$ fb
& $1.42^{+6.9\%}_{-5.1\%}$ fb
\\
\hline
\end{tabular}
\caption{Fiducial cross section of
  $pp\rightarrow\textrm{HZ}{}\rightarrow\left(b\bar{b}\right)\left(e^{+}e^{-}\right)$
  at 13 TeV with leptonic and $b$-jet cuts. The uncertainty band
  refers to the scale variation described in the text. Numerical
  errors for each prediction are beyond the quoted digits.
   }
\label{tab:fiducial-xs}
\end{table}
The fiducial cross sections in this phase-space volume at different
levels of our simulations, are reported in
Tab.~\ref{tab:fiducial-xs}. We also present results with an additional
cut on the Z boson transverse momentum, $p_{t,\rm{Z}} > 150$ GeV,
which we refer to as high-$p_{t,\rm{Z}}$ region.

We first discuss the results without $gg\!\to\!\rm{HZ}{}$ contribution,
over the full range of Z boson transverse momentum reported in
the first line of the Tab.~\ref{tab:fiducial-xs}. The
\noun{HZJ-MiNLO}{} cross section
is about 8\% smaller than the full NNLO calculation from \texttt{MCFM-8.0}{}. The
difference is properly accounted for by reweighting the event sample
and the cross section of \noun{HZ-NNLOPS(LHEF)}{} and \texttt{MCFM-8.0}{} are equal to each
other within the numerical accuracy (which is at the level of the last
quoted digit). The scale uncertainty from the NLO result is reduced
from about 7\% to below 1\% at the NNLO level.
The inclusion of the $\mathcal{O}(\alpha_{\scriptscriptstyle \mathrm{s}}^2)$ $gg\!\to\!\rm{HZ}{}$ channel, reported
in the second line of the table, results in further increase of the
total cross section by about 10\%. In this case, the scale uncertainty
is dominated by the new contribution, which is described only at
leading order, and increases the scale uncertainty to the level of
2-3\%. This larger scale uncertainty is somehow welcome, as a scale
uncertainty below the percent level is unlikely to reflect the true
perturbative uncertainty.  This uncertainty will be reduced by an NLO
treatment of the $gg\!\to\!\rm{HZ}{}$ contribution.\footnote{Note that the small
  difference in the $gg\!\to\!\rm{HZ}{}$ contribution in \texttt{MCFM-8.0}{} and \noun{HZ-NNLOPS(LHEF)}{}
  is due to using \noun{LHAPDF}{} or \noun{POWHEG}{} routines to perform the
  running of the coupling from $M_{\rm{Z}}$ to the central scale choice
  $M_{\rm{H}}+M_{\rm{Z}}$.}

We now discuss the impact of the parton shower on these cross sections. As
is well known, in the presence of fiducial cuts that constrain the jet
activity, as is in the case at hand, there can be a sizeable
difference between a pure fixed-order computation and results after
applying a parton shower. This is illustrated in the last two columns
of the table.
The parton shower allows for extra QCD radiation off coloured partons
which can move the $b$-jets outside the fiducial phase-space volume,
thereby reducing the recorded cross section.  The impact of parton
shower is similar in both instances, with and without the $gg\!\to\!\rm{HZ}{}$ contribution, and amounts to about $10$\% reduction of the cross section in the fiducial region,
while the impact is milder, $5-7$\%, in the high-$p_{\rm t,Z}$ region.

If we now examine the results with an additional $p_{t, \rm{Z}}$ cut,
reported in the last two lines of the table, we observe a reduction of
the cross section by a factor of about 5 and in general behaviours
similar to the ones described above. One point to note is that the
impact of the $gg\!\to\!\rm{HZ}{}$ contribution is larger in this phase space
region, which implies also larger scale uncertainties.

To further illustrate the effect of the $gg\!\to\!\rm{HZ}{}$ channel, we present
in Fig.~\ref{fig:gghz-plots} the differential distributions of the
invariant mass of the \textrm{HZ}{} system and transverse momentum of the
$b$-jet pair system ($p_{\rm t, b\bar b}$) associated to the
reconstruction of the Higgs boson momentum. To define the
\textrm{HZ}{}-invariant mass we use the Monte Carlo truth, while $p_{\rm t,b
  \bar b}$ is obtained by clustering events with the flavour-$k_t$
algorithm with $R=0.4$ and by summing the transverse momenta of the
two $b$-jets.  If more than two $b$-flavoured jets are found, one
selects the pair whose invariant mass is closest to the Higgs
invariant mass.\footnote{We do not distinguish between $b$ and $\bar
  b$ jets.}
We show results from the \noun{HZNNLOPS}{} simulation before
and after the inclusion of the $gg\!\to\!\rm{HZ}{}$ contribution (blue and red
respectively) together with the fixed-order NNLO prediction including
the $gg\!\to\!\rm{HZ}{}$ contribution (green).
For the invariant mass distribution, the large impact of the $gg\!\to\!\rm{HZ}{}$
contribution above the top threshold is evident. Similarly, the
transverse momentum distribution is mostly affected by the $gg\!\to\!\rm{HZ}{}$
contribution in the region between 150 and 200 GeV. In both cases the
impact of $gg\!\to\!\rm{HZ}{}$ remains large up to high scales.

The invariant mass distribution (left panel of
Fig.~\ref{fig:gghz-plots}) features an almost uniform shift between
the fixed-order predictions (green) and the ones including parton
shower evolution (red).  As discussed in the previous paragraph,
parton shower leads to a decrease in the fiducial cross section mainly
due to the $b$-jet cuts. However, the \textrm{HZ}{} invariant mass is not
strictly correlated with the Higgs kinematics (and hence with the
$b$-jets from its decay).
As a consequence, we observe a moderate and constant difference
between \texttt{MCFM-8.0}{} and \noun{HZNNLOPS}{}.

On the other hand, the effect of the parton shower on the transverse
momentum of the reconstructed Higgs boson (right plot of
Fig.~\ref{fig:gghz-plots}) is quite different. In this case, one can
see that the parton shower has a sizeable effect for
$p_{t,b\bar{b}{}}\gtrsim 120$ GeV, and that it smears the distribution in
a non-uniform way. At larger values of $p_{t,b\bar{b}{}}\gtrsim 300$ GeV,
the effect of the partons shower becomes much more modest.
A more detailed discussion is presented in
App.~\ref{App:gghz}.

\begin{figure}
  \centering
  \includegraphics[page=1,width=205pt]{graphics/pheno-new/mhz-pythia-r04-ptz-tot-fid.pdf}
  \hspace*{0.25cm}
  \includegraphics[page=1,width=205pt]{graphics/pheno-new/ptbb-pythia-r04-ptz-tot-fid.pdf}
  \caption{The differential distributions of the invariant mass of the
    \textrm{HZ}{} system (left panel) and the transverse momentum of the Higgs
    boson reconstructed from $b$-jets (right panel). The lower panel
    illustrate ratio of full results (NNLO as well as NNLOPS) to the
    NNLOPS results without $gg\!\to\!\rm{HZ}{}$ contribution.}
  \label{fig:gghz-plots}
\end{figure}

As a next step, we want to study the quality of the Higgs boson
reconstruction.  In Fig.~\ref{fig:higgs-reconstr},
\begin{figure}
  \centering
  \includegraphics[page=1,width=205pt]{graphics/pheno-new/ptbb-reconstruction-r04-ptz-tot-fid-gghz0.pdf}
  \hspace*{0.25cm}
  \includegraphics[page=1,width=205pt]{graphics/pheno-new/ptbb-reconstruction-r04-ptz-tot-fid-gghz1.pdf}
  \includegraphics[page=1,width=205pt]{graphics/pheno-new/ptbb-reconstruction-r07-ptz-tot-fid-gghz0.pdf}
  \hspace*{0.25cm}
  \includegraphics[page=1,width=205pt]{graphics/pheno-new/ptbb-reconstruction-r07-ptz-tot-fid-gghz1.pdf}
  \caption{The differential distributions of the transverse momentum
    of the Higgs boson. \texttt{MC-truth} label refers to the actual
    Higgs boson momentum as passed from the event generator, other
    lines represent the reconstruction of the Higgs boson momentum
    using the two $b$-jets with invariant mass closest to
    $M_{\rm{H}}{}$. The upper two plots show results for jets clustered
    with $R=0.4$, the lower plots $R=0.7$. Left plots do not include
    the $gg\!\to\!\rm{HZ}{}$ contribution, while the right plots do. }
    \label{fig:higgs-reconstr}
\end{figure}
we present a comparison of the
transverse momentum distribution of the true Higgs boson, obtained
using the momentum passed from the event generator before it splits to
$b$-quarks and before any radiation off the $b$-quarks (labelled as
\texttt{MC-truth}), and the $b\bar{b}{}$-jet system, reconstructed with a
flavour-$k_t$ algorithm as described above. We compare two sets of plots
obtained with a jet radius of $R=0.4$ (upper plots) or $R=0.7$ (lower
plots).\footnote{We note that the fiducial cuts are applied on jets of
  $R=0.7$ in this case.} The plots show a comparison of the
    fixed-order results (green curve)\footnote{Note that here only the
      $b$-quarks from the Higgs decay are considered flavoured.}, the
    \noun{HZNNLOPS}{} after parton shower evolution (red) and the MC-truth
    prediction obtained with \noun{HZNNLOPS}{} after parton shower
    (blue). The baseline in the ratio plot is taken to be the latter.
In the left panels the $gg\!\to\!\rm{HZ}{}$ contribution is not included, whereas
its effect is included in right panels.

We start by examining the results without $gg\!\to\!\rm{HZ}{}$ contribution (left
hand panels). We note that both the fixed-order (green) and the
\noun{HZNNLOPS}{} after parton shower (red) differ from the
$\texttt{MC-truth}$ result (blue). At low transverse momenta, this
difference becomes smaller when a larger jet-radius is considered
(left bottom panel), which suggests that the dominant reason for the
difference is out-of-jet radiation from the $b\bar b$-final state. At larger
transverse momenta the difference with respect to the
$\texttt{MC-truth}$ is instead smaller at smaller jet-radius (top left
panel), which points to the fact that in this region the difference is
mainly due to radiation from the initial state.
We also notice that in the intermediate transverse momentum region the
fixed-order and \noun{HZNNLOPS}{} show sizeable differences for small jet
radius, while these differences are more moderate when using a larger
$R$. This can be easily understood from the fact that the
observable with larger $R$ is more inclusive and hence fixed-order and
parton shower results are in better agreement.

We now move to discuss the plots including the $gg\!\to\!\rm{HZ}{}$ effects.
First, we note that the red and green bands in the top right panel if
Fig.~\ref{fig:higgs-reconstr} are identical to the bands shown in the
right panel of Fig.~\ref{fig:gghz-plots}.
As expected when the radius becomes bigger (bottom right panel) the
fixed-order (green) and parton shower results (red) move closer to
each other, again because the observables become more inclusive.
We also note that the uncertainty bands are now larger compared to the
results without $gg\!\to\!\rm{HZ}{}$ contribution. This was already observed for
the fiducial cross section and is due to the leading order
description of the $gg\!\to\!\rm{HZ}{}$ contribution.

We now show the distribution of the transverse momentum of the
$b\bar{b}{}$-jet system in the fiducial volume with and without the
additional cut $p_{\rm t,Z} > 150$ GeV. The relevant plots are shown
in Fig.~\ref{fig:ptbb-jets-zoom}. First of all we note that the
difference between treating the ${\textrm{H}}\!\to\!b\bar{b}{}$ decay at NLO with respect to
LO is very small, which leads to the conclusion that a parton shower
equipped with Matrix Element corrections to the ${\textrm{H}}\!\to\!b\bar{b}{}$ branching
provides a very good estimation of the higher-order corrections. We
also notice a Sudakov shoulder in the fixed-order prediction in the
right panel of Fig.~\ref{fig:ptbb-jets-zoom} at $p_{{\rm t},b\bar b} =
150$ GeV. This feature has already been observed in Figs.~(6) and (12)
of ref.~\cite{Caola:2017xuq} and is due to the fact that the presence
of the $p_{\rm t,Z} > 150$ GeV cut makes the differential spectrum
sensitive to soft gluon emission close to the cut. As expected, a parton
shower captures parts of the resummation effects and therefore the
shoulder is not present in the \noun{NNLOPS}{} predictions.
\begin{figure}
  \centering
  \includegraphics[page=1,width=205pt]{graphics/pheno-new/ptbb-pythia-r04-ptz-tot-fid-gghz1-zoom.pdf}
  \hspace*{0.25cm}
  \includegraphics[page=1,width=205pt]{graphics/pheno-new/ptbb-pythia-r04-ptz-high-fid-gghz1-zoom.pdf}
  \caption{The differential distribution of the transverse momentum of
    the $b\bar{b}{}$-jet system without (left) and with (right) the cut
    $p_{\rm t,Z} > 150$ GeV. Results include the $gg\!\to\!\rm{HZ}{}$
    contribution.
  }
  \label{fig:ptbb-jets-zoom}
\end{figure}

One of the most important variables when reconstructing a resonance is
the invariant mass of its decay products, therefore we will focus on
it in the following, in the boosted high-$p_{\rm t,Z}$ region.
At LO in the decay the $M_{b\bar b}$ distribution is an extremely narrow
Breit-Wigner function, and receives sizeable corrections away from the
peak only at higher-orders.
We start by examining how well \PYTHIA{8} can describe the decay of
the Higgs boson by comparing two calculations that include LO or NLO
decay in the matrix element. When the matrix element is computed at LO
only, \PYTHIA{8} performs the shower using with hardness scale given in
the Les Houches event file.
This comparison is shown in Fig.~\ref{fig:mbb-LOvsNLO} without (left
plot) and with $gg\!\to\!\rm{HZ}$ (right plot). We compare \noun{HZNNLOPS}{} with LO
treatment of the Higgs decay (purple), \noun{HZNNLOPS}{} with NLO
corrections to the ${\textrm{H}}\!\to\!b\bar{b}{}$ decay (red) and \noun{HZJ-MiNLO}{} predictions,
with NLO decay (green).
We see that the two \noun{HZNNLOPS}{} predictions are compatible with each
other all the way down to relatively low $M_{b\bar b}$ masses. We note that
the scale uncertainty is very small, of the order of 2-5\%, when no
$gg\!\to\!\rm{HZ}{}$ contribution is included. This uncertainty increases when
$gg\!\to\!\rm{HZ}{}$ events are included, since these events sit at $M_{b\bar b} = M_{\rm{H}}$
before showering. The small scale variation band is not indicative of
the true uncertainty on this distribution and is related to the fact
that \noun{HZJ-MiNLO}{} results have been reweighted to NNLO results. In
fact, the pure \noun{HZJ-MiNLO}{} predictions, even without $gg\!\to\!\rm{HZ}{}$, have a
larger uncertainty. We also note that this uncertainty is also somehow
underestimated as the band does not cover the \noun{HZNNLOPS}{}
results. This is related to the well known fact that, in a plain
\noun{POWHEG}{} simulation, the scale is varied at the level of the $\bar B$
function, which is by definition inclusive over radiation, whereas the
$M_{b\bar b}$ spectrum is sensitive to radiation.

\begin{figure}
  \centering
  \includegraphics[page=1,width=205pt]{graphics/pheno-new/mbb-pythia-r04-ptz-high-fid-dec-comparison-gghz0.pdf}
  \includegraphics[page=1,width=205pt]{graphics/pheno-new/mbb-pythia-r04-ptz-high-fid-dec-comparison-gghz1.pdf}
  \caption{The differential distributions of the invariant mass of the
    $b\bar{b}{}$-system used for reconstruction of the Higgs
    boson. Comparison of \noun{HZNNLOPS}{} results with LO and NLO decay
    matrix elements, excluding $gg\!\to\!\rm{HZ}{}$ channel (left panels) and
    with $gg\!\to\!\rm{HZ}{}$ (right panels).}
  \label{fig:mbb-LOvsNLO}
\end{figure}

In Fig.~\ref{fig:mbb-plots} we now compare fixed-order predictions
(green) and our best prediction \noun{HZNNLOPS}{} with NLO corrections to
the ${\textrm{H}}\!\to\!b\bar{b}{}$ decay (red).
\begin{figure}
  \centering
  \includegraphics[page=1,width=205pt]{graphics/pheno-new/mbb-pythia-r04-ptz-high-fid-gghz0.pdf}
  \hspace*{0.25cm}
  \includegraphics[page=1,width=205pt]{graphics/pheno-new/mbb-pythia-r04-ptz-high-fid-gghz1.pdf}
  \includegraphics[page=1,width=205pt]{graphics/pheno-new/mbb-pythia-r07-ptz-high-fid-gghz0.pdf}
  \hspace*{0.25cm}
  \includegraphics[page=1,width=205pt]{graphics/pheno-new/mbb-pythia-r07-ptz-high-fid-gghz1.pdf}
  \caption{The differential distributions of the invariant mass of the
    $b\bar{b}{}$-system used for reconstruction of the Higgs boson. We
    present the results obtained with jet clustering with $R=0.4$ (top) and $R=0.7$ (bottom)
    excluding $gg\!\to\!\rm{HZ}{}$ channel (left panels) and with $gg\!\to\!\rm{HZ}{}$ (right panels).}
  \label{fig:mbb-plots}
\end{figure}
In the plots of Fig.~\ref{fig:mbb-plots} we
show predictions obtained with $b$-jets clustered with $R=0.4$ (top
panels) and $R=0.7$ (bottom panels).
We point out that in order to populate the region to the left of the
peak ($M_{b\bar b} < M_{\rm{H}}$) there must be a radiation off the
$b$-quarks produced in the Higgs decay. On the contrary, the region on
the right hand side of the peak ($M_{b\bar b} > M_{\rm{H}}$) is filled
only when additional radiation, off the partons from the production stage,
is clustered with the Higgs decay products.

In Fig.~\ref{fig:mbb-plots} we notice a sizeable enhancement in the
$M_{b\bar b}$ distribution to the left of the Higgs peak. This enhancement
was already observed in refs.~\cite{Ferrera:2017zex,Caola:2017xuq} and
is even more dramatic in this case.
If we compare our left plots to the Figs. (4) and (11) of
ref.~\cite{Caola:2017xuq} we observe a larger K-factor. However there
are a number of differences. First, the results of
ref.~\cite{Caola:2017xuq} are obtained with $R=0.5$. Second, our
\texttt{MCFM-8.0}{} predictions are obtained using massive $b$-quarks, while the
NNLO-approx calculation shown in ref.~\cite{Caola:2017xuq} is obtained
with massless $b$-quarks. Furthermore, the two computations use
different fiducial cuts and in~\cite{Caola:2017xuq} the process
\textrm{HW}{}$^-$ is considered, rather than \textrm{HZ}{}. Last, our plots show
\noun{HZNNLOPS}{} results rather than \noun{HZJ-MiNLO}{} ones, and from
Fig.~\ref{fig:mbb-LOvsNLO} this amounts to a further increase of the
ratio by 10\% (25\%) without (with) $gg\!\to\!\rm{HZ}{}$.
Similar considerations apply when comparing to Fig. (2) of
ref.~\cite{Ferrera:2017zex}.

By looking at the plots on the right of Fig.~\ref{fig:mbb-plots}, one
can observe an even more pronounced enhancement of the \noun{HZNNLOPS}{}
over the \texttt{MCFM-8.0}{} K-factor when the $gg\!\to\!\rm{HZ}{}$ contribution is
included. This is again due to the fact that the $gg\!\to\!\rm{HZ}{}$ term in the
fixed-order calculation is only in the $M_{b\bar b}= M_{\rm{H}}$ bin, while this
contribution is spread to other bins by the parton shower.
A second observation is that when a large jet radius is considered
(bottom row), more radiation is clustered in the $b$-jets. As a
consequence, the distribution vanishes faster away from $M_{\rm{H}}$. This
effect is stronger when a parton shower is included and causes the
K-factor to be smaller to the left of the peak and to even become
close to one at very low mass.

\begin{figure}
  \centering
  \includegraphics[page=1,width=205pt]{graphics/pheno-new/mbb-pythia-all-r04-ptz-high-fid-gghz0-fine-binning.pdf}
  \includegraphics[page=1,width=205pt]{graphics/pheno-new/mbb-pythia-all-r04-ptz-high-fid-gghz1-fine-binning.pdf}
  \caption{The differential distributions of the invariant mass of the
    $b\bar{b}{}$-system used for reconstruction of the Higgs
    boson. Comparison of \noun{HZNNLOPS}{} results with NLO decay matrix
    elements at parton-, hadron-level and with multi-parton
    interactions (MPI), excluding $gg\!\to\!\rm{HZ}{}$ channel (left panels) and
    with $gg\!\to\!\rm{HZ}{}$ (right panels).}
  \label{fig:mbb-NLO-pythia-comparison}
\end{figure}

Finally, in Fig.~\ref{fig:mbb-NLO-pythia-comparison} we present the
$M_{b\bar b}$ distribution obtained with \noun{HZNNLOPS}{} at various stages after
the \PYTHIA{8} parton showering, namely at parton-level and after
hadronisation, with and without multi-parton interactions (MPI). We
notice that hadronisation smears the distribution close to the
peak. This is the reason for the dip at $M_{b\bar b}=M_{\rm{H}}$ in the first inset.
Away from the peak $M_{b\bar b}=M_{\rm{H}}$, we observe that hadronisation effects
become more important at low invariant masses, while, as expected,
they become negligible at large $M_{b\bar b}$. Since color-connections tend
to reduce spatial distances between partons during hadronisation, more
radiation is clustered within the $b$-jets. This is the reason why the
small $M_{b\bar b}$ regions are underpopulated with respect to predictions at
parton level. This effect is similar with or without MPI.  On the
contrary, we can see a substantial change when considering MPI in the
region $M_{b\bar b} > M_{\rm{H}}$. Since MPI provide more radiation activity, many
additional hadrons can be clustered within the $b$-jets, thereby
increasing the invariant mass of the $b\bar{b}{}$-system and causing
migration of events from the region $M_{b\bar b}\approxM_{\rm{H}}$ to larger
invariant masses.

\section{Conclusions}
\label{sec:conclu}
In this paper we have implemented a consistent matching of NNLO
accurate predictions for \textrm{HZ}{} production to parton shower, including
the subsequent decay of the Z boson into pair of leptons and the NLO
decay of the Higgs boson into pair of $b$-quarks.  The \noun{HZNNLOPS}{}
generator we obtained allows for a fully-exclusive simulation of the
\textrm{HZ}{} production in a hadronic collision maintaining the advantages of
the NNLO fixed-order calculation and supplying it with resummation
effects as provided by the matching to a parton shower.

In order to obtain this accuracy, we have extended the existing
\noun{HZJ-MiNLO}{} implementation to include the NLO corrections to the
${\textrm{H}}\!\to\!b\bar{b}{}$ decay. The NNLO+PS matching procedure requires a reweighting
of the \noun{HZJ-MiNLO}{} events that is fully-differential in the \textrm{HZ}{} Born
kinematics. By using properties of the matrix elements at hand, we have
parametrised the latter using variables that allow to express the
differential cross section in terms of a finite set of functions,
thereby simplifying considerably the multi-differential reweighting
procedure. We have also included in the simulation the loop-induced
$gg\!\to\!\rm{HZ}{}$ channel that enters at $\mathcal{O}(\alpha_{\scriptscriptstyle \mathrm{s}}^2)$, as it
constitutes a sizeable part of the total cross section, and can give
rise to substantial distortions of kinematic distributions.

In Sec.~\ref{sec:pheno}, we have considered a setup similar to the one
used in searches for the Higgs decay into $b$-quarks. We find that
scale uncertainties are substantially reduced when the NNLO
corrections are included. Moreover we notice that the cross section in
the fiducial region is reduced by about 5-10\% during parton shower
evolution as a consequence of requiring two $b$-jets satisfying the
fiducial cuts. This correction brings the final result outside the
NNLO uncertainty band. This highlights the limitation of using the
scale uncertainty as an estimate of the true theoretical error
associated to missing higher orders, in particular when more exclusive
fiducial cuts are applied. A NNLO+PS simulation allows to capture some
of these higher order effects, albeit with limited logarithmic
accuracy.

As already noted in the literature, we also find that the $gg\!\to\!\rm{HZ}{}$
channel has a significant impact, especially for the $M_{\sss \HZ} > 2m_t$
part of the spectrum. Moreover its presence has a very strong impact
on the size of the scale uncertainty band. Since this contribution
only enters at $\mathcal{O}(\alpha_{\scriptscriptstyle \mathrm{s}}^2)$ level, in order to reduce this
uncertainty one needs to include higher-order corrections to the
discussed channel, which are currently unknown.

We also notice differences between distributions of the transverse
momentum of the Higgs boson, computed using Monte Carlo truth, and the
transverse momentum of the $b\bar{b}{}$-system identified and used for
the reconstruction of the Higgs boson momentum. We point out that,
especially when a large jet-radius is used, the amount of radiation
clustered into $b$-jets leads to a harder $p_T$-spectrum than the one
of the true Higgs boson.

Despite the consistency of the procedure we used in our simulation, we
have obtained a large K-factor in the $M_{b\bar b}{} < M_{\rm{H}}$ region of the
distribution of the invariant mass of the $b\bar{b}{}$-system. We also
point out that the scale uncertainty of the \noun{NNLOPS}{} prediction in
this part of the spectrum is underestimated, due to known properties
of the algorithm used by \noun{POWHEG}{} to generate real radiation. These
two issues will certainly require further studies. Including NNLO
corrections to the ${\textrm{H}}\!\to\!b\bar{b}{}$ decay and matching them to parton showers
would also be desirable, as well as trying to incorporate NLO
electroweak effects as obtained in
ref.~\cite{Granata:2017iod}. We leave this to future work.

\section*{Acknowledgments}
We would like to thank Gavin Salam for providing a private implementation
of the flavour-$k_t$ algorithm \cite{Banfi:2006hf}. We also are
grateful to Fabrizio Caola and Gionata Luisoni for providing some
results of ref.~\cite{Caola:2017xuq} for cross-checks and to John
Campbell, Giancarlo Ferrera, Keith Hamilton, and Paolo Nason for
useful discussions.
This work was supported in part by ERC Consolidator Grant HICCUP
(No.\ 614577). The work of ER is supported by a Maria
Sk\l{}odowska-Curie Individual Fellowship of the European Commission's
Horizon 2020 Programme under contract number 659147
PrecisionTools4LHC.
WA and WB thank CERN for their hospitality while part of this work was
done.

redpajama_set_name: RedPajamaArXiv