11institutetext: Max Planck Institute for Physics, Föhringer Ring 6, 80805
München, Germany

# NNLO event generation for $pp\to Zh\to\ell^{+}\ell^{-}b\bar{b}$ production
in the SM effective field theory

Ulrich Haisch 1 Darren J. Scott 1 Marius Wiesemann 1 Giulia Zanderighi 1
and Silvia Zanoli<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We consider associated $Zh$ production with $Z\to\ell^{+}\ell^{-}$ and $h\to
b\bar{b}$ decays in hadronic collisions. In the framework of the Standard
Model effective field theory (SMEFT) we calculate the QCD corrections to this
process and achieve next-to-next-to-leading order plus parton shower
(NNLO$+$PS) accuracy using the MiNNLOPS method. This precision is obtained for
a subset of six SMEFT operators, including the corrections from effective
Yukawa- and chromomagnetic dipole-type interactions. Missing higher-order QCD
effects associated with the considered dimension-six operators are estimated
to have a relative numerical impact of less than a percent on the total rate
once existing experimental limits on the relevant Wilson coefficients are
taken into account. We provide a dedicated Monte Carlo (MC) code that
evaluates the NNLO SMEFT corrections on-the-fly in the event generation. This
MC generator is used to study the numerical impact of NNLO$+$PS corrections on
the kinematic distributions in $pp\to Zh\to\ell^{+}\ell^{-}b\bar{b}$
production employing simple SMEFT benchmark scenarios. We identify the
invariant mass $m_{b\bar{b}}$ of the two $b$-tagged jets as well as the three-
invariant jet mass $m_{b\bar{b}j}$ as particularly interesting observables to
study SMEFT effects. These distributions receive contributions that change
both their normalisation and shape with the latter modifications depending on
the exact jet definition. To our knowledge SMEFT effects of this type have so
far not been discussed in the literature. The presented MC generator can also
serve as a starting point to obtain NNLO$+$PS accuracy for a suitable enlarged
set of effective operators in the future.

##  1 Introduction

The Standard Model (SM) predicts that a Higgs boson with a mass of $125\,{\rm
GeV}$ decays with a branching ratio close to 60% to a pair of bottom quarks.
The most sensitive production channels for detecting $h\to b\bar{b}$ decays at
the LHC are despite their notably smaller cross sections with respect to
gluon-gluon fusion Higgs production, the associated production of a Higgs
boson with a $W$ or $Z$ boson ($Vh$), where the leptonic decay of the vector
boson enables a clean selection. The $h\to b\bar{b}$ decay mode has been
observed by both the ATLAS and the CMS collaboration in LHC Run II
ATLAS:2018kot ; CMS:2018nsn , and these measurements constrain the $h\to
b\bar{b}$ signal strength in $Vh$ production
$\big{(}\mu_{b\bar{b}}^{Vh}\big{)}$ to be SM-like within about 25%, at the
level of one standard deviation. With ATLAS and CMS being able to collect new
data very soon in LHC Run III and with the high-luminosity upgrade (HL-LHC) on
the horizon, the precision of the $\mu_{b\bar{b}}^{Vh}$ measurements will
improve significantly with an ultimate projected HL-LHC accuracy of $15\%$
($5\%$) in the case of $Wh$ ($Zh)$ production ATLAS:2018jlh ; CMS:2018qgz .

Future LHC measurements of in particular $pp\to Zh\to\ell^{+}\ell^{-}b\bar{b}$
production will thus allow to place strong constraints on non-standard Higgs
interactions. Assuming that any potential beyond-the-SM (BSM) contributions
arises only from particles with masses much heavier than the electroweak (EW)
scale, the SM effective field theory (SMEFT) Buchmuller:1985jz ;
Grzadkowski:2010es ; Brivio:2017vri  offers a largely model-independent,
systematically improvable quantum-field theoretical framework for analyses of
BSM effects in collider processes that can be used to probe Higgs production
and decay. In fact, radiative corrections in the SMEFT to both $Vh$ production
Mimasu:2015nqa ; Degrande:2016dqg  and the $h\to f\bar{f}$ decays
Gauld:2015lmb ; Gauld:2016kuu ; Cullen:2019nnr ; Cullen:2020zof  have been
calculated. The recent paper Bizon:2021rww  has furthermore studied QCD
corrections to $pp\to Zh\to\ell^{+}\ell^{-}b\bar{b}$ in the anomalous-coupling
framework. The existing studies for $Vh$ production have focused on the subset
of higher-dimensional interactions that modify the couplings of the Higgs to
EW gauge bosons achieving next-to-leading order (NLO) Mimasu:2015nqa ;
Degrande:2016dqg  and next-to-next-to-leading order (NNLO) Bizon:2021rww  in
QCD perturbation theory, respectively, while in the case of $h\to b\bar{b}$
both NLO QCD and NLO EW corrections to the total decay width have been
calculated for the full set of relevant dimension-six SMEFT operators
Gauld:2015lmb ; Gauld:2016kuu ; Cullen:2019nnr .

The ultimate goal of a SMEFT calculation of the $pp\to
Zh\to\ell^{+}\ell^{-}b\bar{b}$ process consists in taking into account all
relevant contributions of dimension-six operators to both $pp\to Zh$ and $h\to
b\bar{b}$, and to combine the production and the decay processes into a Monte
Carlo (MC) event generator that consistently includes higher-order QCD and EW
corrections at fixed order as well as parton-shower (PS) effects. Within the
SM such computations have reached NNLO$+$PS precision Astill:2018ivh ;
Alioli:2019qzz ; Bizon:2019tfo ; Zanoli:2021iyp , which means that they
include NNLO QCD corrections to $pp\to Zh\to\ell^{+}\ell^{-}b\bar{b}$
production as well as the matching to PS MC generators. A dedicated PS MC code
including both NLO QCD and EW corrections within the SM also exists
Granata:2017iod .

The intent of this work is to achieve NNLO$+$PS accuracy in the SMEFT
including all dimension-six operator insertions that contribute to the
subprocesses $pp\to Zh$ and $h\to b\bar{b}$ directly in QCD. This requires to
calculate QCD corrections for a subset of six SMEFT operators, including
effects associated to effective Yukawa- and chromomagnetic dipole-type
interactions. The inclusion of these contributions in the NNLO+PS event
generation yields an accurate description of the SMEFT effects in differential
predictions for the full $pp\to Zh\to\ell^{+}\ell^{-}b\bar{b}$ reaction. Since
the distribution of the invariant mass $m_{b\bar{b}}$ of the two $b$-tagged
jets ($b$-jets) receives contributions that change both its normalisation and
shape it turns out to provide a particularly useful observable in the context
of the SMEFT. The latter modifications depend on the $b$-jet definition, a
feature that has not been discussed in the SMEFT literature as far as we know.
Similar observations can be made for the three-jet invariant mass
$m_{b\bar{b}j}$. As a by-product of our computations we are also able to
extend the calculation of the inclusive $h\to b\bar{b}$ decay rate in the
SMEFT Gauld:2016kuu  to the next order in QCD for the case that bottom quarks
are treated as massless. On the other hand, non-trivial EW corrections arising
in the SMEFT, such as a modified $hZZ$ coupling, are not taken into account in
our analysis. Since QCD and EW SMEFT corrections to $pp\to
Zh\to\ell^{+}\ell^{-}b\bar{b}$ approximately factorise, in the sense that the
EW operators contributing to production play only a subleading role in the
decay and vice versa, obtaining NNLO$+$PS accuracy for a suitable enlarged set
of effective operators should be rather straightforward. While beyond the
scope of this publication, we plan to return to this problem in the future.

This manuscript is structured as follows: in Section 2 we specify the subset
of dimension-six operators of the full SMEFT Lagrangian that are relevant in
the context of this article, including a discussion of the normalisation
chosen for the individual effective interactions. Section 3 contains a brief
description of the basic ingredients of the SMEFT calculations for $pp\to Zh$
and $h\to b\bar{b}$ and their combination and implementation in our NNLO$+$PS
event generator. The impact of the SMEFT corrections on kinematic
distributions in $pp\to Zh\to\ell^{+}\ell^{-}b\bar{b}$ production at NNLO$+$PS
is presented in Section 4 by using simple benchmark scenarios for the Wilson
coefficients. We conclude and present an outlook in Section 5. The lenghty
analytic expressions for the squared matrix elements that are relevant for our
work are relegated to Appendix A, while Appendix B contains numerical
estimates of higher-order QCD corrections associated to the subset of the
SMEFT operators that are considered in this paper. The discussed corrections
have been neglected in our phenomenological study because they all turn out to
contribute less than a percent once existing experimental limits on the
relevant Wilson coefficients are taken into account.

##  2 Preliminaries

In this article we consider the following set of dimension-six operators

$\displaystyle Q_{H\Box}$ $\displaystyle=(H^{\dagger}H)\hskip
1.42262pt\Box\hskip 1.42262pt(H^{\dagger}H)\,,$ $\displaystyle Q_{HD}$
$\displaystyle=(H^{\dagger}D_{\mu}H)^{\ast}\hskip
1.42262pt(H^{\dagger}D^{\mu}H)\,,$ $\displaystyle Q_{bH}$
$\displaystyle=y_{b}\hskip 0.7113pt(H^{\dagger}H)\hskip
1.42262pt\bar{q}_{L}\hskip 0.7113ptb_{R}\hskip 0.7113ptH\,,$ $\displaystyle
Q_{bG}$ $\displaystyle=\frac{g_{s}^{3}}{(4\pi)^{2}}\hskip 1.42262pty_{b}\hskip
1.42262pt\bar{q}_{L}\sigma_{\mu\nu}T^{a}\hskip 0.7113ptb_{R}\hskip
0.7113ptH\hskip 0.7113ptG^{a,\hskip 0.7113pt\mu\nu}\,,$ (1) $\displaystyle
Q_{HG}$ $\displaystyle=\frac{g_{s}^{2}}{(4\pi)^{2}}\hskip
1.42262pt(H^{\dagger}H)\hskip 1.42262ptG_{\mu\nu}^{a}\hskip
0.7113ptG^{a,\hskip 0.7113pt\mu\nu}\,,$ $\displaystyle Q_{3G}$
$\displaystyle=\frac{g_{s}^{3}}{(4\pi)^{2}}\hskip 1.42262ptf^{abc}\hskip
0.7113ptG_{\mu}^{a,\hskip 0.7113pt\nu}G_{\nu}^{b,\hskip
0.7113pt\sigma}G_{\sigma}^{c,\hskip 0.7113pt\mu}\,,$

which appear in the full SMEFT Lagrangian

${\cal L}_{\rm SMEFT}\supset\sum_{i}\frac{C_{i}}{\Lambda^{2}}\,Q_{i}\,.$ (2)

Here $\Box=\partial_{\mu}\hskip 0.7113pt\partial^{\mu}$,
$\sigma_{\mu\nu}=i/2\hskip
1.42262pt(\gamma_{\mu}\gamma_{\nu}-\gamma_{\nu}\gamma_{\mu})$ with
$\gamma_{\mu}$ the usual Dirac matrices, $H$ denotes the SM Higgs doublet,
$q_{L}$ is the left-handed third-generation quark doublet, $b_{R}$ is the
right-handed bottom-quark singlet, while $g_{s}=\sqrt{4\pi\alpha_{s}}$ and
$G_{\mu\nu}^{a}$ denote the coupling constant and the field strength tensor of
QCD, respectively. The definition of the covariant derivative is
$D_{\mu}=\partial_{\mu}-i\hskip 0.7113ptg_{s}\hskip 0.7113ptG_{\mu}^{a}\hskip
0.7113ptT^{a}$ with $T^{a}$ being the $SU(3)$ generators and $f^{abc}$ denote
the fully antisymmetric QCD structure constants. The bottom-quark Yukawa
coupling is defined as $y_{b}=\sqrt{2}\hskip 0.7113pt\bar{m}_{b}/v$, with the
$\overline{\rm MS}$ bottom-quark mass $\bar{m}_{b}$ and the Higgs vacuum
expectation value (VEV) $v$, while $\Lambda$ denotes the new-physics mass
scale that suppresses the dimension-six operators $Q_{i}$ entering (2) and
$C_{i}$ are the corresponding Wilson coefficients. Notice finally that in the
case of $Q_{bH}$ and $Q_{bG}$ the sum over the hermitian conjugate in (2) is
understood.

The normalisations of the dimension-six operators introduced in (2) deserve
some additional comments. First, the two mixed-chirality operators $Q_{bH}$
and $Q_{bG}$ include a factor of $y_{b}$ which serves as an order parameter
and explicitly appears in a broad class of ultraviolet (UV) completions that
match onto the set of operators in (2). See for example the discussions in
Giudice:2007fh ; Elias-Miro:2013mua . Second, the factors of $g_{s}$ and
$1/(4\pi)^{2}$ that arise in the definition of $Q_{bG}$, $Q_{HG}$ and $Q_{3G}$
guarantee that the associated Wilson coefficients $C_{bG}$, $C_{HG}$ and
$C_{3G}$ are expected to be of ${\cal O}(1)$ in all weakly-coupled UV-complete
extensions of the SM with new degrees of freedom and masses in the ballpark of
$\Lambda$.111In the case of the operator $Q_{bG}$ the corresponding Wilson
coefficient can also be of ${\cal O}(g^{2}/g_{s}^{2})$ with $g$ a weak
coupling which implies that $C_{bG}$ is parametrically smaller than ${\cal
O}(1)$. This happens when $Q_{bG}$ is generated by a EW and not a strong loop
in the UV theory. Notice that in the operator basis that has been employed in
the NLO QCD calculation of the inclusive $h\to b\bar{b}$ decay rate in the
SMEFT Gauld:2016kuu  a different normalisation is chosen for the operators
$Q_{bH}$, $Q_{bG}$ and $Q_{HG}$. As a result, in this normalisation the
corresponding Wilson coefficients are expected to be of size $C_{bH}={\cal
O}(y_{b})={\cal O}\big{(}10^{-2}\big{)}$, $C_{bG}={\cal O}\big{(}y_{b}\hskip
0.7113pt\alpha_{s}/(4\pi)\big{)}={\cal O}\big{(}10^{-4}\big{)}$ and
$C_{HG}={\cal O}\big{(}\alpha_{s}/(4\pi)\big{)}={\cal O}\big{(}10^{-2}\big{)}$
in weakly-coupled BSM theories and not $C_{i}={\cal O}(1)$ as in the operator
basis (2). Our normalisation therefore has the merit that the suppression
factors $y_{b}$ and $\alpha_{s}/(4\pi)$ appear explicitly as order parameters
which allows for a more explicit power counting in our SMEFT calculation of
QCD corrections to the $pp\to Zh\to\ell^{+}\ell^{-}b\bar{b}$ process.

##  3 Calculation in a nutshell

In this section, we describe the different ingredients of the calculation of
the QCD corrections to the fully differential decay rate of $h\to b\bar{b}$
and the $Zh$ production cross section in the SMEFT. Throughout this work, we
use the five-flavour scheme and thus treat the bottom quark as massless both
in the matrix elements and the phase-space integrals. The bottom-quark Yukawa
coupling is however taken to be non-zero. The explicit expressions for the
non-trivial $h\to b\bar{b}$ squared matrix elements can be found in Appendix
A. Moreover, we assume minimal-flavour violation DAmbrosio:2002vsn  and set
the Cabibbo-Kobayashi-Maskawa matrix element $V_{tb}$ to unity. After having
discussed the anatomy of the SMEFT corrections to both $pp\to Zh$ production
and the $h\to b\bar{b}$ decay, we briefly mention the employed NNLO$+$PS
methods and explain how we apply them to the event generation of the $pp\to
Zh\to\ell^{+}\ell^{-}b\bar{b}$ process including SMEFT effects.

###  3.1 Factorisable contributions to the $h\to b\bar{b}$ decay

Since the operators $Q_{H\Box}$, $Q_{HD}$ and $Q_{bH}$ do not contain a gluon
the associated SMEFT contributions factorise to all orders in $\alpha_{s}$. As
a result, the matrix elements proportional to the Wilson coefficients
$C_{H\Box}$, $C_{HD}$ and $C_{bH}$ can be obtained from the massless NNLO
calculation of the fully differential $h\to b\bar{b}$ decay rate within the SM
Anastasiou:2011qx ; DelDuca:2015zqa ; Caola:2019pfz  by the following simple
replacement:

$y_{b}^{2}\to y_{b}^{2}\hskip 2.13394pt\big{(}1+2\hskip 0.35565ptc_{\rm
fac}\big{)}\,,$ (3)

with

$c_{\rm fac}=c_{\rm kin}-c_{bH}\,,\qquad c_{\rm
kin}=\frac{v^{2}}{\Lambda^{2}}\left[C_{H\Box}-\frac{C_{HD}}{4}\right]\,,\qquad
c_{bH}=\frac{v^{2}}{\Lambda^{2}}\hskip 1.42262pt{\rm
Re}\left(C_{bH}\right)\,.$ (4)

Notice that the term $c_{\rm kin}$ arises from the canonical normalisation of
the Higgs kinetic term in the presence of $Q_{H\Box}$ and $Q_{HD}$. The
squared Higgs VEV entering (4) is related to the Fermi constant $G_{F}$
extracted from muon decay via $v^{2}=1/\big{(}\sqrt{2}\hskip
0.7113ptG_{F}\big{)}$. EW tree-level corrections appearing in the SMEFT that
modify the relation between the VEV and the Fermi constant $\big{(}$see for
example Brivio:2017vri ; Gauld:2015lmb $\big{)}$ are very small and hence
neglected in our analysis. In practice, our implementation relies on the
expressions for the squared matrix elements provided in DelDuca:2015zqa  which
also have been used in the publications Astill:2018ivh ; Bizon:2019tfo ;
Zanoli:2021iyp  to obtain NNLO$+$PS predictions for $pp\to
Zh\to\ell^{+}\ell^{-}b\bar{b}$ production within the SM.

Given the simplicity of the replacement (3) it is straightforward to obtain an
analytic result for the factorisable corrections to the inclusive $h\to
b\bar{b}$ decay rate due to $Q_{H\Box}$, $Q_{HD}$ and $Q_{bH}$ up to NNLO in
QCD. The leading-order (LO) expression for the partial Higgs boson decay width
to massless bottom-quark pairs within the SM is given by

$\Gamma(h\to b\bar{b})_{\rm SM}^{\rm LO}=\frac{3\hskip 0.7113pty_{b}^{2}\hskip
0.7113ptm_{h}}{16\pi}\,.$ (5)

Using this expression the corresponding NNLO result takes the form
Gorishnii:1983cu ; Gorishnii:1990zu ; Gorishnii:1991zr

$\Gamma(h\to b\bar{b})_{{\rm SM}}^{\rm NNLO}=\big{(}1+\Delta\big{)}\hskip
1.42262pt\Gamma(h\to b\bar{b})_{\rm SM}^{\rm
LO}\,,\qquad\Delta=\frac{\alpha_{s}}{\pi}\hskip
1.42262pt5.67+\left(\frac{\alpha_{s}}{\pi}\right)^{2}\hskip 1.42262pt29.15\,.$
(6)

Here both the bottom-quark Yukawa coupling $y_{b}$ and the strong coupling
constant $\alpha_{s}$ are understood to be renormalised at the scale $m_{h}$
and we have assumed five active quark flavours to obtain the numerical results
for the expansion coefficients appearing in $\Delta$. By means of the
replacement (3) one then finds

$\Gamma(h\to b\bar{b})_{{\rm SMEFT}}^{{\rm NNLO},{\rm fac}}=\big{(}1+2\hskip
0.35565ptc_{\rm fac}\big{)}\hskip 1.42262pt\Gamma(h\to b\bar{b})_{{\rm
SM}}^{\rm NNLO}\hskip 0.7113pt\,,$ (7)

where the expressions for $c_{\rm fac}$ and $\Gamma(h\to b\bar{b})_{{\rm
SM}}^{\rm NNLO}$ have already been given in (4) and (6), respectively. Since
in the limit of massless bottom quarks the inclusive $h\to b\bar{b}$ decay
rate is known to ${\cal O}(\alpha_{s}^{4})$ Chetyrkin:1996sr ; Baikov:2005rw ;
Herzog:2017dtz  an extension of (7) to ${\rm N}^{4}{\rm LO}$ in QCD would be
possible. Such higher-order QCD corrections are included in eHDECAY
Contino:2013kra ; Contino:2014aaa , which allows to calculate the factorisable
SMEFT corrections to $h\to b\bar{b}$ in the so-called strongly-interacting
light Higgs or SILH Giudice:2007fh  basis of dimension-six operators. Since
the aim of this paper is to consistently achieve NNLO$+$PS accuracy for the
fully differential decay rate of $h\to b\bar{b}$ in the SMEFT, the result in
(7) will however be sufficient for our purposes.

###  3.2 Non-factorisable contributions to the $h\to b\bar{b}$ decay

Insertions of the operator $Q_{bG}$ lead to non-factorisable QCD corrections
to the fully differential $h\to b\bar{b}$ decay rate. Examples of the
respective Feynman diagrams are shown in Figure 1. The relevant squared matrix
elements are found by interfering the SMEFT contributions with the
corresponding SM amplitudes. The explicit expressions for the resulting
squared matrix elements can be found in Appendix A. The leading contribution
arises at ${\cal O}\big{(}y_{b}^{2}\hskip 0.7113pt\alpha_{s}^{2}\hskip
0.7113ptC_{bG}\big{)}$ from the interference of the $h\to b\bar{b}g$ amplitude
in the SMEFT (upper left diagram) and the corresponding SM graphs. Notice that
this real contribution is IR finite and that the corresponding one-loop
contribution to $h\to b\bar{b}$ is identical to zero since it only involves
scaleless integrals. In the operator basis (2) the first non-factorisable QCD
corrections involving the operator $Q_{bG}$ therefore appear at NNLO in QCD.
In terms of the inclusive LO SM decay rate (5), we find by integrating (23)
over the three-particle phase space the following compact expression:

$\Gamma(h\to b\bar{b})_{\rm SMEFT}^{{\rm NNLO},{\rm non}}=\Delta_{\rm
non}\hskip 1.42262ptc_{bG}\hskip 1.42262pt\Gamma(h\to b\bar{b})_{\rm SM}^{\rm
LO}\,,\qquad\Delta_{\rm non}=\left(\frac{\alpha_{s}}{\pi}\right)^{2}\hskip
1.42262pt\frac{m_{h}^{2}}{3\hskip 0.7113ptv^{2}}\,,$ (8)

with

$c_{bG}=\frac{v^{2}}{\Lambda^{2}}\hskip 1.42262pt{\rm
Re}\left(C_{bG}\right)\,.$ (9)

We emphasise that when taking the difference in the normalisation of the
operator $Q_{bG}$ into account, the result (8) agrees with the expression
derived in Gauld:2016kuu  after taking the limit $m_{b}\to 0$.

Figure 1:  Examples of SMEFT contributions to the fully differential $h\to
b\bar{b}$ decay rate involving an insertion of the operator $Q_{bG}$ (black
square). The upper left (right) diagram represents a tree-level (one-loop)
contribution to the $h\to b\bar{b}g$ decay, the center left (right) diagram
yields a tree-level contribution to the $h\to b\bar{b}q\bar{q}$ $\big{(}$$h\to
b\bar{b}gg$$\big{)}$ process, while the lower diagrams contribute to the $h\to
b\bar{b}$ amplitude at the two-loop level. Notice that the quark flavours in
the center left $h\to b\bar{b}q\bar{q}$ diagram can be $q=u,d,s,c,b$ and that
effective five-point $hb\bar{b}gg$ vertices also contribute in the case of the
$h\to b\bar{b}gg$ transition. See text for further details.

In Section 2 we have argued that the Wilson coefficient of $Q_{bG}$ is
expected to be of ${\cal O}(1)$ in a wide class of UV-complete theories if the
operator is normalised as in (2). While there are strong bounds on the
imaginary part of $C_{bG}$ from the electric dipole moment of the neutron
Haisch:2021hcg , the real part entering (9) is at present only very weakly
constrained by experiments Hayreter:2013kba ; Bramante:2014hua ;
Haisch:2021hcg . In fact, values of $c_{bG}={\cal O}(100)$ are compatible with
all existing low- and high-energy data (see Section 4 for details), and
therefore from a purely phenomenological point of view it is possible that
$\alpha_{s}/\pi\hskip 1.42262ptc_{bG}={\cal O}(1)$ or even larger. In such a
case the NNLO correction in (8) is numerically of the size of a NLO
correction. To deal with this possibility we decided to include in our work
also all contributions of ${\cal O}\big{(}y_{b}^{2}\hskip
0.7113pt\alpha_{s}^{3}\hskip 0.7113ptC_{bG}\big{)}$. As shown in Figure 1
there are three types of such contributions. There are virtual (upper right
diagram) or real corrections to the Born-level process $h\to b\bar{b}g$, the
latter of which lead to either $h\to b\bar{b}q\bar{q}$ (center left diagram)
or $h\to b\bar{b}gg$ (center right diagram). Notice that the four-quark final
state involves the quark flavours $q=u,d,s,c,b$ and that final-state
configurations with two bottom quarks and two gluons receive contributions
from both effective $hb\bar{b}g$ and $hb\bar{b}gg$ vertices. Also the two-loop
virtual corrections to the $h\to b\bar{b}$ decay (lower diagrams) contribute
at ${\cal O}\big{(}y_{b}^{2}\hskip 0.7113pt\alpha_{s}^{3}\hskip
0.7113ptC_{bG}\big{)}$. In Section 4 we will show that the inclusion of the
N3LO corrections associated to $Q_{bG}$ are phenomenologically relevant, which
provides a further rational to incorporate them in our NNLO$+$PS calculation
of $pp\to Zh\to\ell^{+}\ell^{-}b\bar{b}$ production.

We now turn our attention to the matrix elements involving the insertion of an
operator $Q_{HG}$ or $Q_{3G}$. In the former case one finds that in the limit
of a massless bottom quark the operator $Q_{HG}$ does not contribute to the
$h\to b\bar{b}$ decay. This is related to the fact that for $m_{b}=0$ an
insertion of $Q_{HG}$ does not lead to amplitudes such as $h\to b\bar{b}$
where the final-state bottom quarks have mixed chiralities. This feature is
well-known from the SM calculation of the inclusive $h\to b\bar{b}$ rate where
amplitudes involving a $h\to gg$ subdiagram first contribute at NNLO in QCD
Chetyrkin:1995pd ; Larin:1995sqc  and only if non-zero quark masses are
considered. In fact, the latter results receive Sudakov-like double logarithms
of the form $\alpha_{s}^{2}\ln^{2}(m_{b}^{2}/m_{h}^{2})$ from both real and
virtual contributions. The same feature is present in the SMEFT calculations
of $h\to f\bar{f}$ processes Gauld:2016kuu ; Cullen:2019nnr ; Cullen:2020zof .
While bottom-quark mass effects in the fully differential SM $h\to b\bar{b}$
decay rate have been studied at NNLO Bernreuther:2018ynm ; Behring:2019oci ;
Behring:2020uzq ; Somogyi:2020mmk , it is not known how to resum the double-
logarithmic contributions in the case of the $h\to b\bar{b}$ decay.222In the
case of the inclusive $h\to\gamma\gamma$ and $h\to gg$ decays it has been
shown how to resum the leading double-logarithmic corrections Akhoury:2001mz ;
Liu:2017vkm ; Liu:2019oav ; Wang:2019mym . Likewise, it is also unknown how to
correctly treat such terms in any of the existing NNLO$+$PS procedures. We
therefore neglect all corrections associated to the insertion of $Q_{HG}$,
which is formally correct in the limit of massless bottom quarks. We show in
Appendix B that the corrections to the inclusive $\Gamma(h\to b\bar{b})$ decay
rate proportional to $C_{HG}$ and at leading power in $m_{b}$ can be bounded
in a model-independent fashion to the level of a few permille. Neglecting such
effects in the calculation of the fully differential $\Gamma(h\to b\bar{b})$
decay rate can therefore also be assumed to be an excellent approximation.

Employing the operator basis introduced in (2) insertions of $Q_{3G}$ induce
tree-level corrections to $h\to b\bar{b}gg$ and one-loop corrections to $h\to
b\bar{b}g$. After interfering these two types of contributions with the
relevant SM amplitudes the resulting SMEFT corrections are proportional to
${\cal O}\big{(}y_{b}^{2}\hskip 0.7113pt\alpha_{s}^{3}\hskip
0.7113ptC_{3G}\big{)}$ and hence formally of N3LO. Since our goal is it to
only attain NNLO$+$PS precision we neglect corrections due to $Q_{3G}$. In
Appendix B we present an estimate of this type of N3LO corrections that
suggests that given the existing bounds on the Wilson coefficient $C_{3G}$,
SMEFT contributions of ${\cal O}\big{(}y_{b}^{2}\hskip
0.7113pt\alpha_{s}^{3}\hskip 0.7113ptC_{3G}\big{)}$ can indeed only have a
very minor numerical impact on the $h\to b\bar{b}$ decay distributions.

###  3.3 Contributions to $pp\to Zh$ production

Let us first recall that the only effective dimension-six EW interactions that
we are taking into account in our analysis are related to the three operators
$Q_{H\Box}$, $Q_{HD}$ and $Q_{bH}$ introduced in (2) and that otherwise we
only consider effective interactions that directly induce QCD corrections. In
particular, dimension-six SMEFT operators that lead to non-trivial
modifications of the $hZZ$ vertex are not considered in our work.

With this simplification the dominant corrections to $Zh$ production in the
SMEFT are associated to $Q_{H\Box}$ and $Q_{HD}$. These operators already
provide a contribution at Born level which can be obtained by the shift

$g_{hZZ}^{2}\to g_{hZZ}^{2}\hskip 2.13394pt\big{(}1+2\hskip 0.35565ptc_{\rm
kin}\big{)}\,,$ (10)

where $g_{hZZ}=2\hskip 0.35565ptm_{Z}^{2}/v$ denotes the SM coupling between a
Higgs and two $Z$ bosons with $m_{Z}$ the $Z$-boson mass. The coefficient
$c_{\rm kin}$ has been defined in (4). Applying the rescaling (10) to the
inclusive $Zh$ production cross section in the SM one obtains the following
formula

$\sigma(pp\to Zh)_{\rm SMEFT}^{\rm NNLO}=\big{(}1+2\hskip 0.35565ptc_{\rm
kin}\big{)}\hskip 1.42262pt\sigma(pp\to Zh)_{\rm SM}^{\rm NNLO}\,,$ (11)

for the SMEFT corrections to the cross section up to NNLO. An analogous
expression also holds at the differential level.

Let us also discuss the role of operators other than $Q_{H\Box}$ and $Q_{HD}$
in associated Higgs production with a vector boson. For massless bottom quarks
neither $Q_{bH}$ nor $Q_{bG}$ furnishes a non-zero contribution to $pp\to Zh$
to all orders in the strong coupling constant $\alpha_{s}$ due to the mixed-
chirality nature of the two operators. Insertions of $Q_{HG}$ lead to a tree-
level and an one-loop correction to $q\bar{q}\to Zhg$ and $q\bar{q}\to Zh$,
respectively. After interfering these amplitudes with their SM counterparts
one obtains a contribution that is of ${\cal O}(\alpha_{s}^{2}\hskip
1.42262ptC_{HG})$ with respect to the LO $Zh$ production cross section within
the SM. While these corrections are therefore formally needed to achieve
NNLO$+$PS accuracy, we will show in Appendix B that the contributions
proportional to $C_{HG}$ cannot exceed the level of a few permille. This
renders them irrelevant for all practical purposes and we hence neglect them.
We note that by using the results of Brein:2011vx  their inclusion would be
quite straightforward, but given their subleading impact we refrain from doing
so. Finally, the operator $Q_{3G}$ gives rise to tree-level (one-loop)
corrections to $q\bar{q}\to Zhgg$ ($q\bar{q}\to Zhg$). Both contributions are
of ${\cal O}(\alpha_{s}^{3}\hskip 1.42262ptC_{3G})$ and therefore we do not
take them into account in our analysis. In Appendix B we nevertheless estimate
their potential size and show them to be negligible for our purposes.

###  3.4 NNLO$+$PS calculation and MC implementation

In the following, we discuss how the fixed-order QCD results including the set
of dimension-six operators in (2) are implemented into the NNLO$+$PS accurate
$h\to b\bar{b}$ generator developed in Bizon:2019tfo . We then describe how
the $h\to b\bar{b}$ decay events can be combined with Higgs production events
from any Higgs production mode, including SMEFT effects consistently. While in
this paper we consider the case of associated $Zh$ production, the
modifications to the formulae required for other processes are straightforward
as long as one works in the narrow width approximation for the Higgs
propagator.

We start our discussion by briefly recalling how NNLO$+$PS accurate results
for the $h\to b\bar{b}$ decay can be achieved in the SM. The simplest way to
obtain NNLO$+$PS precision is to compute decay events by applying the so-
called MiNLO′ method Hamilton:2012np ; Hamilton:2012rf . The weights are then
rescaled to the NNLO accurate inclusive $h\to b\bar{b}$ decay rate
$\Gamma(h\to b\bar{b})_{\rm SM}^{\rm NNLO}$ by multiplying the event weights
with the ratio $\Gamma(h\to b\bar{b})_{\rm SM}^{\rm NNLO}/\Gamma(h\to
b\bar{b})_{\rm SM}^{{\rm MiNLO}^{\prime}}$, where the denominator is directly
computed from the MiNLO′ events by summing their weights. As a variant, it is
also possible to restrict the NNLO correction to a certain phase space region.
In particular, one can use a rescaling factor that tends to one in the region
where the decay events involve hard radiation, while requiring the integral of
the rescaled events to reproduce exactly $\Gamma(h\to b\bar{b})_{\rm SM}^{\rm
NNLO}$. In this work, we follow such an approach and specifically employ the
procedure described in Section 2.2 of the article Bizon:2019tfo .

As discussed in Sections 3.1 and 3.2, the SMEFT corrections to $h\to b\bar{b}$
split into factorisable and non-factorisable corrections. The factorisable
SMEFT effects can be taken into account by using an analogous reweighting to
(7) and applying it to the NNLO$+$PS accurate $h\to b\bar{b}$ events, i.e. by
computing $\big{(}1+2\hskip 0.35565ptc_{\rm fac}\big{)}\hskip
0.7113pt\Gamma(h\to b\bar{b})_{\rm SM}^{\rm NNLO}/\Gamma(h\to b\bar{b})_{\rm
SM}^{{\rm MiNLO}^{\prime}}$ when reweighting the MiNLO′ events. As far as non-
factorisable corrections are concerned, they are separately IR finite, as
pointed out already before. In particular, the leading contribution including
$Q_{bG}$ is of ${\cal O}(\alpha_{s}^{2})$ and contributes to the squared
matrix element of $h\to b\bar{b}g$, cf. (23), which features no divergence
when the gluon is unresolved and will be referred to as ${\rm R}$ hereafter.
Hence, the integration over the three-body phase space is IR finite without
the need to apply the MiNLO′ method, at variance with the NNLO$+$PS
calculation of the factorisable part. At ${\cal O}(\alpha_{s}^{3})$ the one-
loop virtual and real corrections to the $h\to b\bar{b}g$ process are
included, just as in a standard POWHEG NLO$+$PS calculation Nason:2004rx ;
Frixione:2007vw . In the following, we will refer to these contributions as
${\rm RV}$ and ${\rm RR}$, respectively. Since their sum is IR finite, these
corrections yield a contribution to the inclusive $h\to b\bar{b}$ decay rate
when integrated over the three- and four-body phase space, respectively, and
added together. To deal with the soft and collinear singularities of the real
contributions and to cancel the IR poles of the one-loop virtual corrections,
cf. (26), we exploit the general implementation of the Frixione-Kunszt-Signer
(FKS) subtraction Frixione:1995ms ; Frixione:1997np  within the POWHEG-BOX
framework Alioli:2010xd . To this end, the full POWHEG-BOX machinery is used
that automatically builds the soft and collinear counterterms and remnants,
and also checks the behaviour in the soft and collinear limits of the real
squared matrix elements against their soft and collinear approximations.
Finally, also the two-loop virtual corrections to $h\to b\bar{b}$ including
$Q_{bG}$ contribute to the non-factorisable SMEFT effects at ${\cal
O}(\alpha_{s}^{3})$. These corrections, called ${\rm VV}$ in what follows, are
IR finite by themselves, cf. (40), and after integrating them over the two-
particle phase space we find the compact expression

$\Gamma(h\to b\bar{b})_{\rm SMEFT}^{{\rm non},{\rm
VV}}=\frac{\alpha_{s}}{\pi}\hskip 0.35565pt\frac{49}{12}\hskip
1.42262pt\Delta_{\rm non}\hskip 1.42262ptc_{bG}\hskip 1.42262pt\Gamma(h\to
b\bar{b})_{\rm SM}^{\rm LO}\,,$ (12)

with $\Delta_{\rm non}$ and $c_{bG}$ defined as in (8) and (9), respectively.
The result in (12) can simply be added to the factorisable SMEFT contribution
to the $h\to b\bar{b}$ decay.

To explain how we combine production and decay events in the SMEFT to obtain
$pp\to Zh\to\ell^{+}\ell^{-}b\bar{b}$ events, it is useful to first recall how
such a combination can be achieved within the SM. In this case, the Higgs
boson in each $Zh$ production event is replaced by the possible decay products
(i.e. $b\bar{b}g$, $b\bar{b}q\bar{q}$ or $b\bar{b}gg$) taken from a given
decay event. Working in the narrow width approximation, the weight $w_{\rm
full}^{\rm SM}$ of the full SM event is then calculated as

$w_{\rm full}^{\rm SM}=\frac{w_{\rm prod}^{\rm SM}\,w_{\rm dec}^{\rm
SM}}{\Gamma_{h}^{\rm SM}}\,,$ (13)

where $w_{\rm prod}^{\rm SM}$ denotes the weight of the production event
obtained at NNLO$+$PS using the MiNNLOPS $pp\to Zh$ generator presented in
Zanoli:2021iyp , while $w_{\rm dec}^{\rm SM}$ represents the weight of the
decay event computed at NNLO$+$PS precision employing the MiNLO′ method and
procedure detailed in the work Bizon:2019tfo . The total decay width
$\Gamma_{h}^{\rm SM}$ of the $125\,{\rm GeV}$ SM Higgs boson serves as an
input. The Les Houches event (LHE) file also stores, for each event, the value
of the hardest radiation allowed by the shower (i.e. the scale scalup). We
record in the combined LHE file the value of scalupprod for the production
process and, once the event is passed to the PS, we recompute on-the-fly the
value of scalupdec for the specific decay kinematics. We then generate
emissions of the Higgs decay in all the available phase space and after the
shower is complete we check whether the hardness of the splittings is below
the veto scale scalupdec. If this is not the case we attempt to shower the
event again until the condition is met. Additional details of this procedure
are given in Section 3.2 of the paper Bizon:2019tfo .

Figure 2:  Illustration of the interference contribution that leads to the
SMEFT correction proportional to $c_{HG}$ in the partial decay width $h\to gg$
as given in (15). The left and right diagram represent the tree-level SMEFT
with an $Q_{HG}$ insertion (black square) and the one-loop SM contribution
involving a top-quark loop, respectively. Light-quark loops also appear in the
SM, but since these corrections vanish identically for $m_{q}=0$, they turn
out to be numerically insignificant and hence are ignored in our study.

From the above discussion one can deduce that NNLO$+$PS accurate differential
cross sections for the full process $pp\to Zh\to\ell^{+}\ell^{-}b\bar{b}$ in
the SMEFT can be obtained by combining production and decay as follows

$\begin{split}d\sigma_{\text{NNLO$+$PS}}&=\big{(}1+2\hskip 0.35565ptc_{\rm
kin}\big{)}^{2}\,\Bigg{\\{}\left[1-2\hskip 0.35565ptc_{bH}+\frac{\Gamma(h\to
b\bar{b})_{\rm SMEFT}^{{\rm non},{\rm VV}}}{\Gamma(h\to b\bar{b})_{{\rm
SM}}^{\rm NNLO}}\right]\hskip 0.7113ptd\sigma^{\rm
SM}_{\text{NNLO$+$PS}}\\\\[5.69054pt] &\hskip 78.24507pt\,+d\sigma^{\rm
non,{\rm R}+{\rm RV}+{\rm RR}}_{\text{NNLO$+$PS}}\hskip
0.7113pt\Bigg{\\}}\;\frac{\Gamma_{h}^{\rm SM}}{\Gamma_{h}^{\rm
SMEFT}}\,,\end{split}$ (14)

where we have factorised the term $\big{(}1+2\hskip 0.35565ptc_{\rm
kin}\big{)}^{2}$ that arises from the canonical normalisation of the Higgs
kinetic term. The results for $\Gamma(h\to b\bar{b})_{{\rm SM}}^{\rm NNLO}$
and $\Gamma(h\to b\bar{b})_{\rm SMEFT}^{{\rm non},{\rm VV}}$ can be found in
(6) and (12), $d\sigma^{\rm SM}_{\text{NNLO$+$PS}}$ represents the NNLO$+$PS
accurate differential cross section obtained interfacing the MiNNLOPS
calculation of the $pp\to Zh$ generator Zanoli:2021iyp  with the reweighted
MiNLO′ calculation of the $h\to b\bar{b}$ generator Bizon:2019tfo , while
$d\sigma^{\rm non,{\rm R}+{\rm RV}+{\rm RR}}_{\text{NNLO$+$PS}}$ includes the
non-factorisable corrections ${\rm R}$, ${\rm RV}$ and ${\rm RR}$ that we
compute as discussed above. Notice that both $\Gamma(h\to b\bar{b})_{\rm
SMEFT}^{{\rm non},{\rm VV}}$ and $d\sigma^{\rm non,{\rm R}+{\rm RV}+{\rm
RR}}_{\text{NNLO$+$PS}}$ depend linearly on the Wilson coefficient $c_{bG}$.
We stress that the advantage of expressing (14) in this form lies in the fact
that arbitrary combinations of Wilson coefficients can be obtained without
recalculating any of the individual cross sections. Thus, variations of the
Wilson coefficients can be obtained a posteriori.

The last factor in (14) takes into account that the total decay width of the
Higgs boson that appears in (13) is modified by SMEFT effects. In our
implementation we employ the following result

$\begin{split}\Gamma_{h}^{\rm SMEFT}&=\big{(}1+2\hskip 0.35565ptc_{\rm
kin}\big{)}\Bigg{[}\,\Gamma_{h}^{\rm SM}-\big{(}2\hskip 0.7113pt\Delta\hskip
0.7113ptc_{bH}-K_{bG}\hskip 1.42262pt\Delta_{\rm non}\hskip
0.7113ptc_{bG}\big{)}\hskip 1.42262pt\Gamma(h\to b\bar{b})_{\rm SM}^{\rm
LO}\\\\[5.69054pt] &\hskip 71.13188pt+6\hskip 0.7113ptK_{HG}\hskip
0.7113ptc_{HG}\hskip 0.7113pt\Gamma(h\to gg)_{\rm SM}^{\rm
LO}\,\Bigg{]}\,.\end{split}$ (15)

The LO expression for the partial decay width for $h\to b\bar{b}$ can be found
in (5) and the corrections proportional to $\Gamma(h\to b\bar{b})_{\rm
SM}^{\rm LO}$ have been included in (15) in the approximation that treats the
bottom quark as strictly massless. The relevant correction factors $\Delta$
and $\Delta_{\rm non}$ have been defined in (6) and (8), respectively. The
multiplicative factor $K_{bG}$ encodes the QCD corrections up to N3LO related
to the $Q_{bG}$ contribution to the partial decay width of $h\to b\bar{b}$.
The used numerical value of $K_{bG}=1.622$ follows from the semi-analytic
formula given below in (21). The term in $\Gamma_{h}^{\rm SMEFT}$ proportional
to

$c_{HG}=\frac{v^{2}}{\Lambda^{2}}\hskip 1.42262ptC_{HG}\,,$ (16)

encodes the SMEFT corrections to the partial decay width of $h\to gg$
involving a single insertion of the operator $Q_{HG}$. The relevant diagrams
are shown in Figure 2. These terms start at ${\cal O}(\alpha_{s}^{2}\hskip
0.7113ptC_{HG})$ since

$\Gamma(h\to gg)^{\rm LO}_{\rm SM}=\frac{\alpha_{s}^{2}\hskip
0.7113ptm_{h}^{3}}{72\hskip 0.7113pt\pi^{3}\hskip 0.7113ptv^{2}}\,.$ (17)

This result holds in the limit of an infinite top-quark mass, which is a
sufficiently accurate approximation for our purpose. To include higher-order
effects proportional to $C_{HG}$ in an approximate way we have incorporated
the factor $K_{HG}=1.844$ in (15). This $K$-factor includes QCD corrections up
to N4LO Herzog:2017dtz . Notice finally that (15) contains only QCD
corrections related to the dimension-six SMEFT operators introduced in
(2).333By adding the term $\Delta\Gamma_{h}^{\rm SMEFT}=-\big{(}1+2\hskip
0.35565ptc_{\rm kin}\big{)}\sum_{q=u,d,s,c}\big{(}2\hskip 0.7113pt\Delta\hskip
0.7113ptc_{qH}-K_{qG}\hskip 0.7113pt\Delta_{\rm non}\hskip
0.7113ptc_{qG}\big{)}\hskip 1.42262pt\Gamma(h\to q\bar{q})_{\rm SM}^{\rm LO}$
to (15) with $K_{qG}=1.622$ it is straightforward to include the dominant
SMEFT corrections to the partial decay widths for $h\to q\bar{q}$.

We finally recall from the discussion in Section 3.3 that if one considers
associated $Zh$ production alone, the SMEFT production weight would be
modified by a universal factor leading to

$w_{\rm prod}^{\rm SMEFT}=\left(1+2\hskip 0.35565ptc_{\rm kin}\right)w_{\rm
prod}^{\rm SM}\,.$ (18)

The multiplicative factor $\left(1+2\hskip 0.35565ptc_{\rm kin}\right)$,
however, cancels against the same overall factor appearing in the total decay
width of the Higgs boson when production and decay events are combined. This
cancellation becomes manifest in (14) when (15) is used. If, on the other
hand, one considers production processes other than $pp\to Zh$, where SMEFT
effects give rise to non-universal corrections to the production process, some
of the formulae given in this subsection need to be modified accordingly.

##  4 Phenomenological analysis

In the following we present NNLO$+$PS accurate results for $pp\to
Zh\to\ell^{+}\ell^{-}b\bar{b}$ production at the LHC with a centre-of-mass
energy of $\sqrt{s}=13\,{\rm TeV}$ including the SMEFT effects discussed in
Section 3. All SM input parameters are taken from the most recent PDG review
Zyla:2020zbs . In particular, we use $G_{F}=1.166379\cdot 10^{-5}\,{\rm
GeV}^{-2}$, $m_{Z}=91.1876\,{\rm GeV}$, $\Gamma_{Z}^{\rm SM}=2.4952\,{\rm
GeV}$, $m_{h}=125.09\,{\rm GeV}$ and $\Gamma_{h}^{\rm SM}=4.1\,{\rm MeV}$. The
Higgs VEV is calculated using $v=2^{-1/4}\hskip
0.7113ptG_{F}^{-1/2}=246.22\,{\rm GeV}$, while the value of the
electromagnetic coupling and the weak mixing angle is derived as
$\alpha=\sqrt{2}\hskip 0.7113ptG_{F}\hskip
0.7113ptm_{W}^{2}\sin^{2}\theta_{w}/\pi=1/132.184$ and
$\sin^{2}{\theta_{w}}=1-{m_{W}^{2}}/{m_{Z}^{2}}=0.2230$ employing
$m_{W}=80.379\,{\rm GeV}$. In the case of the bottom-quark mass we use the
$\overline{\rm MS}$ mass $\bar{m}_{b}(\bar{m}_{b})=4.18\,{\rm GeV}$ as input
which leads to $\bar{m}_{b}(m_{h})=2.79\,{\rm GeV}$ and
$y_{b}(m_{h})=\sqrt{2}\hskip 0.7113pt\bar{m}_{b}(m_{h})/v=1.60\cdot 10^{-2}$.
NNPDF31_nnlo_as_0118 parton distribution functions (PDFs) Ball:2017nwa  with
$\alpha_{s}(m_{Z})=0.1180$ giving rise to $\alpha_{s}(m_{h})=0.1127$ are
employed in our MC simulations and events are showered with Pythia 8
Sjostrand:2014zea  utilising the Monash tune Skands:2014pea . Effects from
hadronisation, underlying event modelling or QED effects in the shower are not
included.

We study two SMEFT benchmark scenarios. In both cases we set the Wilson
coefficients $c_{\rm kin}$, $c_{HG}$ and $c_{3G}$ to zero — see (4) and (16).
In the first scenario, we choose $c_{bH}\neq 0$ and $c_{bG}=0$, while in the
second one we take $c_{bH}=0$ and $c_{bG}\neq 0$. Before specifying our
benchmark settings for the two parameters, let us summarise the existing
experimental constraints on $c_{bH}$ and $c_{bG}$. In the former case, we take
the result from a recent global SMEFT fit Ellis:2020unq  that includes 34
dimension-six operators. In the normalisation of (2) and (4), this work
provides the following marginalised 95% confidence level (CL) bound

$c_{bH}\in[-0.13,0.20]\,.$ (19)

In the case of $c_{bG}$ one instead has

$c_{bG}\in[-438,438]\,,$ (20)

at 95% CL, which has been obtained in Haisch:2021hcg  from an analysis of the
transverse momentum ($p_{T}$) spectrum of $Z$-boson production in association
with $b$-jets as measured by ATLAS in LHC Run II ATLAS:2020juj . We add that
the nominal strongest bound on $c_{bG}$ has been derived in Haisch:2021hcg
and relies on high-mass $b$-jet pair production ATLAS data ATLAS:2019fgd .
While this search imposes $c_{bG}\in[-149,149]$ it is not clear to which
extent a SMEFT treatment is trustworthy in this high-mass region, so that we
do not use it here. The present bounds from Higgs physics on the Wilson
coefficient $c_{bG}$ Hayreter:2013kba ; Bramante:2014hua  are weaker than the
constraint given in (20).

###  4.1 Inclusive $h\to b\bar{b}$ decay in the SMEFT

As a first application of the calculations outlined in Section 3 we extend the
results for the inclusive $h\to b\bar{b}$ decay width in the SMEFT presented
in Gauld:2016kuu  to the N3LO level in QCD for massless bottom quarks. In
terms of the LO SM inclusive decay width (5) we find

$\begin{split}\Gamma(h\to b\bar{b})_{\rm
SMEFT}^{\text{N${}^{3}$LO}}&=\Bigg{\\{}\left(1+2\hskip 0.35565ptc_{\rm
fac}\right)\left[1+\frac{\alpha_{s}}{\pi}\hskip
1.42262pt5.67+\left(\frac{\alpha_{s}}{\pi}\right)^{2}\hskip
1.42262pt29.15+\left(\frac{\alpha_{s}}{\pi}\right)^{3}\hskip
0.7113pt41.76\right]\\\\[5.69054pt]
&\phantom{xxx}+\left(\frac{\alpha_{s}}{\pi}\right)^{2}\frac{m_{h}^{2}}{3\hskip
0.7113ptv^{2}}\left[1+\frac{\alpha_{s}}{\pi}\hskip
1.42262pt17.32\right]c_{bG}\Bigg{\\}}\;\Gamma(h\to b\bar{b})_{\rm
SM}^{\text{LO}}\,,\end{split}$ (21)

if the renormalisation scale is identified with the Higgs boson mass by
setting $\mu=m_{h}$. Notice that the ${\cal O}(\alpha_{s}^{3})$ corrections
proportional to $\left(1+2\hskip 0.35565ptc_{\rm fac}\right)$ are known from
the SM calculation of the $h\to b\bar{b}$ decay Chetyrkin:1996sr ;
Baikov:2005rw ; Herzog:2017dtz . The result for the ${\cal O}(\alpha_{s}^{3})$
correction proportional to $c_{bG}$ is instead new and given here for the
first time. Notice that the latter terms enhance the non-factorisable
contribution due to $Q_{bG}$ by around 60%, which provides a clear motivation
to incorporate them in our NNLO$+$PS generator. The corrections associated
with $Q_{HG}$ and $Q_{3G}$ are very small and not included in (21). See the
discussion in Appendix B.

Using (21) one can now study the possible numerical impact of the set of
dimension-six operators introduced in (2). Allowing for instance $c_{bH}$
($c_{bG}$) to vary within its experimentally allowed 95% CL range in eqs. (19)
and (20), while setting all other Wilson coefficients to zero, leads to the
following relative shifts in the inclusive $h\to b\bar{b}$ decay width:

$\frac{\Gamma(h\to b\bar{b})_{\rm SMEFT}^{\text{N${}^{3}$LO}}}{\Gamma(h\to
b\bar{b})_{\rm
SM}^{\text{N${}^{3}$LO}}}-1\in\left\\{\begin{matrix}[-39,26]\%&\text{for}~{}(\ref{eq:cbH95CL})~{}\text{and}~{}c_{\rm
kin}=c_{bG}=0\,,\\\\[5.69054pt]
[-6.3,6.3]\%&\text{for}~{}(\ref{eq:cbG95CL})~{}\text{and}~{}c_{\rm
kin}=c_{bH}=0\,.\end{matrix}\right.$ (22)

Thus, there is a hierarchy between the possible SMEFT effects in the $h\to
b\bar{b}$ decay rate, with the non-factorisable contributions due to $Q_{bG}$
being smaller by a factor of ${\cal O}(5)$ than the factorisable corrections
that are associated to $Q_{bH}$ (as well as $Q_{H\Box}$ and $Q_{HD}$). The
SMEFT QCD corrections in (43) and (47) due to $Q_{HG}$ and $Q_{3G}$ are
smaller by more than a factor of ${\cal O}(20)$ compared to the contributions
from $Q_{bG}$. Missing higher-order QCD effects in (21) can therefore only
have a relative numerical impact of a few permille once existing experimental
limits on the Wilson coefficients of the operators in (2) are taken into
account.

###  4.2 Differential $pp\to Zh\to\ell^{+}\ell^{-}b\bar{b}$ cross section in
the SMEFT

In our differential analysis we select events with two charged leptons
(electrons or muons) to explore the $Zh\to\ell^{+}\ell^{-}b\bar{b}$ signature.
The leptons are required to have a transverse momentum of $p_{T,\ell}>15\,{\rm
GeV}$ and a pseudorapidity of $|\eta_{\ell}|<2.5$. The invariant mass of the
dilepton pair is restricted to $m_{\ell^{+}\ell^{-}}\in[75,105]\,{\rm GeV}$.
The events are furthermore required to have at least two $b$-jets, which are
reconstructed using the anti-$k_{t}$ algorithm Cacciari:2008gp  as implemented
in FastJet Cacciari:2011ma . We impose transverse momentum cuts of
$p_{T,b}>25\,{\rm GeV}$ and a rapidity threshold of $|\eta_{b}|<2.5$ on the
$b$-jets. The definition of potential additional jets use the same thresholds
as those of the $b$-jets. The dominant background processes are $Z+{\rm
jets}$, $t\bar{t}$, single-top and diboson production. The latter three types
of backgrounds can be substantially reduced by requiring large values of
$p_{T,Z}$ Butterworth:2008iy . Hence, to improve the signal-to-background
ratio we impose $p_{T,Z}\in[150,250]\,{\rm GeV}$. Notice that this $p_{T,Z}$
requirement corresponds to the second resolved $p_{T,Z}$ bin as recommended in
the stage 1.2 simplified template cross sections (STXS) framework
Andersen:2016qtm ; Berger:2019wnu ; Amoroso:2020lgh  which is also implemented
in the latest ATLAS LHC Run II measurements of the $pp\to
Zh\to\ell^{+}\ell^{-}b\bar{b}$ process ATLAS:2019yhn ; ATLAS:2020fcp . We will
also comment on how our results are modified if the other two resolved
regions, i.e. $p_{T,Z}\in[75,150]\,{\rm GeV}$ and $p_{T,Z}>250\,{\rm GeV}$,
are considered.

Figure 3:  Invariant mass of the two $b$-jets using the anti-$k_{t}$ algorithm
with radius parameter $R=0.4$. The red histogram in the left (right) panel
corresponds to the prediction for $c_{bH}=0.15$ ($c_{bG}=400$). For comparison
our SM prediction with its scale uncertainty band is shown in black and gray.
All results correspond to proton-proton ($pp$) collisions at
$\sqrt{s}=13\,{\rm TeV}$ and are subject to the fiducial cuts discussed in the
main text. The lower panels depict the ratios between the BSM and the SM
distributions.

The two panels in Figure 3 display our predictions for the $pp\to
Zh\to\ell^{+}\ell^{-}b\bar{b}$ cross section differential in the invariant
mass of the two $b$-jets, employing a jet radius of $R=0.4$ in the
anti-$k_{t}$ clustering. If more than two $b$-jets are present the observable
$m_{b\bar{b}}$ is defined as the invariant mass of the pair of $b$-jets
closest to $m_{h}$. The black curves correspond to our SM NNLO+PS prediction
for the $13\,{\rm TeV}$ LHC with central renormalisation scale $\mu_{R}$ and
factorisation scale $\mu_{F}$ set according to the MiNNLOPS procedure
Monni:2019whf ; Monni:2020nks  and the gray band represents the corresponding
perturbative uncertainties. These uncertainties have been obtained from seven-
point scale variations enforcing the constraint $1/2\leq\mu_{R}/\mu_{F}\leq 2$
and keeping the scale variation in production and decay correlated. The same
way of estimating perturbative uncertainties is applied to all kinematic
distributions that are provided in this section. The red histogram in the left
and right panel of Figure 3 corresponds to the results for $c_{bH}=0.15$ and
$c_{bG}=400$, respectively. These values are within the range allowed by (19)
and (20). All other Wilson coefficients not specified in a given plot are set
to zero.

From the left plot in Figure 3 it is evident that BSM effects in the form of a
non-zero Wilson coefficient $c_{bH}$ just lead to a rescaling of the
$m_{b\bar{b}}$ spectrum. This is expected because $c_{bH}$ is part of the
factorisable corrections in (4) that just rescales all kinematic distributions
by an overall factor. On the other hand, a non-zero Wilson coefficient
$c_{bG}$ is more interesting, since $c_{bG}\neq 0$ alters the shape of the
$m_{b\bar{b}}$ distribution with respect to the SM prediction. This can be
seen in the right panel of Figure 3. In fact, one observes that for the
choices $c_{bG}=400$ and $R=0.4$ the $m_{b\bar{b}}$ spectrum receives relative
corrections of up to $40\%$ for invariant masses $m_{b\bar{b}}\simeq 50\,{\rm
GeV}$. The reason for this somewhat surprising feature is the structure of the
tree-level squared matrix element, given in (23), that modifies the $h\to
b\bar{b}g$ process and constitutes the leading $Q_{bG}$ contribution. The
corresponding Feynman diagram is shown in Figure 1 in the upper row on the
left-hand side. From (23) one observes that the probability for emitting a
gluon is flat in phase space. In contrast, the real emission contribution to
the differential decay rate $h\to b\bar{b}g$ in the SM is divergent when the
radiated gluon becomes unresolved, i.e. soft or collinear to one of the bottom
quarks, and therefore such emissions are favoured. As a result, configurations
where the total invariant mass $m_{b\bar{b}g}=m_{h}$ of the $b\bar{b}g$ system
is shared equally between the three individual partons occur much more
frequently in the former than in the latter case, where the bottom quarks
typically carry most of the energy which leads to an invariant mass
distribution that is strongly peaked at $m_{b\bar{b}}\simeq m_{h}$. We add
that changing the sign of $c_{bH}$ or $c_{bG}$ will also change the sign of
the relative corrections due to the considered SMEFT operators.

Figure 4:  Same as the right panel in Figure 3, but for $R=0.7$ (left) and
$R=1.0$ (right). See the main text for further details.

Notice that in the case of $c_{bG}\neq 0$ the shape of the $m_{b\bar{b}}$
distribution depends on the jet radius $R$ used to identify $b$-jets. To
illustrate this feature we display in Figure 4 two additional spectra assuming
again $c_{bG}=400$, but taking $R=0.7$ and $R=1.0$ instead of the standard
choice $R=0.4$. One observes that the corrections due to $Q_{bG}$ are on
average pushed towards lower values of $m_{b\bar{b}}$ when the jet radius $R$
is increased. We further add in this context that at ${\cal
O}(\alpha_{s}^{3})$ insertions of the operator $Q_{bG}$ lead to a one-loop
contribution to the $h\to b\bar{b}g$ amplitude, tree-level contributions to
the $h\to b\bar{b}q\bar{q}$ and $h\to b\bar{b}gg$ processes and two-loop
effects in the $h\to b\bar{b}$ amplitude (cf. Figure 1). While the two-loop
corrections contribute only at $m_{b\bar{b}}=m_{h}$, the other two types of
contributions are again spread over the phase space. Since the ${\cal
O}(\alpha_{s}^{3})$ non-factorisable corrections due to $Q_{bG}$ are
relatively large in the case of the inclusive $h\to b\bar{b}$ decay width in
(21), including them in the calculation of the differential cross sections for
the full process $pp\to Zh\to\ell^{+}\ell^{-}b\bar{b}$ is necessary if one
wants to describe kinematic distributions such as $m_{b\bar{b}}$ accurately.
We also note that, while the results shown in Figures 3 and 4 have been
obtained for $p_{T,Z}\in[150,250]\,{\rm GeV}$, qualitatively similar
modifications of the $m_{b\bar{b}}$ distribution due to $Q_{bG}$ are found
when the transverse momentum of the $Z$ boson is restricted to the other two
stage 1.2 STXS regions with $p_{T,Z}\in[75,150]\,{\rm GeV}$ and
$p_{T,Z}>250\,{\rm GeV}$. As for the shapes, the observed differences depend
on the exact $b$-jet definition. In fact, they turn out to be more pronounced
for a larger radius $R$.

Figure 5:  Same as Figure 4, but for $m_{b\bar{b}j}$ using $R=0.4$ (left
panel) and $R=0.7$ (right panel). See the main text for further details.

Another distribution that features interesting shape changes in the presence
of a non-zero Wilson coefficient $c_{bG}$ is the invariant mass
$m_{b\bar{b}j}$ of the two $b$-jets and an extra jet. We built this observable
from the set of two $b$-jets and one additional jet whose three-jet invariant
mass lies closest to the Higgs boson mass. Figure 5 shows two $m_{b\bar{b}j}$
spectra for $c_{bG}=400$ with a jet radius of $R=0.4$ and $R=0.7$ in the left
and the right panel, respectively. Also the modifications in the
$m_{b\bar{b}j}$ spectrum due to $c_{bG}\neq 0$ are non-trivial and
$R$-dependent. In fact, the relative SMEFT effects are more pronounced at
lower values of $m_{b\bar{b}j}$ and can reach up to around 40%. This feature
can be qualitatively understood by remembering that the leading $Q_{bg}$
contributions arise from the process $h\to b\bar{b}g$, see (23), while in the
SM the corresponding matrix element is part of the NLO corrections. Relative
to the inclusive $h\to b\bar{b}$ decay width, events with an additional jet
will therefore occur more likely in SMEFT scenarios with $c_{bG}\neq 0$ than
in the SM. Moreover, since the largest corrections in the invariant mass
distributions of the two $b$-jets arise at $m_{b\bar{b}}\simeq 40\,{\rm GeV}$
(cf. the right panel in Figure 3 and the two panels in Figure 4) and we impose
$p_{T,j}>25\,{\rm GeV}$, one expects to see an excess of events at
$m_{b\bar{b}j}\simeq 60\,{\rm GeV}$. Indeed, this is what is observed in the
left panel of Figure 5. It is also clearly visible from the two plots in
Figure 5 that increasing the jet radius from $R=0.4$ to $R=0.7$ will result in
a migration of events to higher $m_{b\bar{b}j}$ values, since a larger jet
radius will collect more radiation, leading on average to a larger three-jet
invariant mass. We finally mention that in the recent ATLAS analysis
ATLAS:2020fcp  of $pp\to Vh$ production with $h\to b\bar{b}$ decay the mass
$m_{b\bar{b}j}$ of the three-jet system is already used as an input to build
the multivariate discriminant in the case of three-jet events (see Table 5 of
that publication). We believe that besides $m_{b\bar{b}}$ the variable
$m_{b\bar{b}j}$ can play an important role in the context of multivariate
discriminants tailored to put constraints on the Wilson coefficient of the
operator $Q_{bG}$.

Let us add that apart from the $m_{b\bar{b}}$ and $m_{b\bar{b}j}$ spectra we
have identified additional kinematic distributions that are sensitive to the
non-factorisable corrections resulting from the operator $Q_{bG}$. For
instance, also the transverse momentum $p_{T,Z}$ of the $Z$ boson is modified
in a non-trivial fashion by $c_{bG}\neq 0$. However, in view of (22) the
effects in $p_{T,Z}$ cannot exceed the percent level, and therefore this
variable taken by itself will have only a rather limited constraining power at
the LHC. Similar statements apply to the transverse momentum $p_{T,b\bar{b}}$
of the $b$-jet pair. Incorporating the latter observables into a multivariate
discriminant may however enhance the overall sensitivity to BSM effects
associated to $c_{bG}\neq 0$. An analysis of this issue is clearly beyond the
scope of this article. Likewise we also do not attempt to derive bounds on the
Wilson coefficients $c_{bH}$ and $c_{bG}$ using existing ATLAS:2018kot ;
CMS:2018nsn ; ATLAS:2019yhn ; ATLAS:2020fcp  or hypothetical ATLAS:2018jlh ;
CMS:2018qgz  differential LHC data, leaving such an exercise for future
research.

##  5 Conclusions

In this article, we have presented novel predictions for the $pp\to
Zh\to\ell^{+}\ell^{-}b\bar{b}$ process within the SMEFT. NNLO QCD corrections
for both $pp\to Zh$ production and $h\to b\bar{b}$ decay have been calculated
for a subset of six dimension-six operators appearing in the full SMEFT
Lagrangian. These fixed-order predictions have been consistently matched to a
PS, allowing for a realistic exclusive description of the process at the level
of hadronic events while retaining NNLO QCD accuracy. All SMEFT corrections
that are associated to the considered operators and can exceed a percent at
the level of the total rate after taking into account the existing
experimental bounds on the relevant Wilson coefficients have been implemented
into a dedicated MC generator.

Using our NNLO$+$PS generator we have performed a phenomenological study of
the impact of SMEFT contributions on several kinematic distributions in $pp\to
Zh\to\ell^{+}\ell^{-}b\bar{b}$ production considering simple benchmark
scenarios. First, we have employed our calculation to extend the SMEFT
computation of the inclusive $h\to b\bar{b}$ decay width Gauld:2016kuu  to
${\cal O}(\alpha_{s}^{3})$ in the five-flavour scheme. Our result for the N3LO
corrections due to the chromomagnetic dipole-type operator $Q_{bG}$ is new.
These higher-order QCD effects enhance the numerical impact of $Q_{bG}$ by
about 60% compared to the ${\cal O}(\alpha_{s}^{2})$ contributions that were
known before. In our phenomenological analysis of differential distributions
in $pp\to Zh\to\ell^{+}\ell^{-}b\bar{b}$ production we have considered
explicitly the $m_{b\bar{b}}$ and $m_{b\bar{b}j}$ spectra and illustrated how
SMEFT effects can change both the normalisation and the shape of these
distributions. While the factorisable corrections associated to the operators
$Q_{H\Box}$, $Q_{HD}$ and $Q_{bH}$ act as an overall rescaling factor on all
kinematic observables, the non-factorisable contributions that arise from the
operator $Q_{bG}$ in (2) lead to non-trivial modifications on the shapes of
certain distributions. By studying the dependence on the jet radius $R$ that
enters the anti-$k_{t}$ jet clustering algorithm, we have shown that these
shape changes depend on the $b$-jet definition, a feature that to our
knowledge has not been discussed in the SMEFT context before. It might be
possible to exploit the observed $R$-dependence of the kinematic distributions
such as $m_{b\bar{b}}$ or $m_{b\bar{b}j}$ to enhance the LHC sensitivity to
the operator $Q_{bG}$ that is at present only very weakly constrained.
Furthermore, combining the invariant masses $m_{b\bar{b}}$ and $m_{b\bar{b}j}$
with other observables such as $p_{T,Z}$ and $p_{T,b\bar{b}}$ in multivariate
discriminants that describe the full kinematics of the selected events is
likely to enhance the sensitivity of the LHC to BSM realisations with a non-
zero Wilson coefficient $c_{bG}$.444At future $e^{+}e^{-}$ colliders it might
be possible to exploit global event-shape variables in the $h\to b\bar{b}$
decay Coloretti:2022jcl  to put constraints on the Wilson coefficient
$c_{bG}$. Since we believe that our MC generator should prove useful for
everyone interested in comparing LHC data to high-precision SMEFT predictions,
we will make the relevant codes to simulate NNLO$+$PS events for the $pp\to
Zh\to\ell^{+}\ell^{-}b\bar{b}$ process in the SMEFT publicly available within
the POWHEG-BOX framework POWHEGBOXRES .

Our analysis of SMEFT effects in the $pp\to Zh\to\ell^{+}\ell^{-}b\bar{b}$
process includes only a limited subset of dimension-six operators. For
instance effective interactions that modify the couplings of the Higgs to EW
gauge bosons are not considered in this work. Since it is known Mimasu:2015nqa
; Degrande:2016dqg ; Bizon:2021rww  that the latter type of SMEFT
contributions can lead to phenomenological relevant effects in $Zh$
production, we plan to apply the discussed NNLO$+$PS method to non-trivial
insertions of purely EW operators as well. Since the leading EW operators
contributing to $Zh$ production play only a minor role in the $h\to b\bar{b}$
decay and vice versa, obtaining NNLO$+$PS accuracy for a suitable enlarged set
of dimension-six operators in the SMEFT seems possible following the
methodology developed here.

###### Acknowledgements.

We are grateful to Rhorry Gauld for helpful discussions and useful comments
that helped us to improve the manuscript. We also thank Marumi Kado for
interesting conversations. The simulations related to this work have used the
Max Planck Computing and Data Facility (MPCDF) in Garching. Silvia Zanoli is
supported by the International Max Planck Research School (IMPRS) on
“Elementary Particle Physics”.

##  Appendix A Squared matrix elements

Below we provide the expressions for the squared matrix elements for the non-
factorisable SMEFT terms of ${\cal O}(y_{b}^{2}\hskip 0.7113ptC_{bG})$ to the
$h\to b\bar{b}$ decay distributions that are included in our NNLO$+$PS
generator. All matrix elements are computed using conventional dimensional
regularisation for both UV and IR singularities in $d=4-2\hskip
0.35565pt\epsilon$ dimensions with $\overline{\rm MS}$ subtraction. The actual
generation and computation of squared matrix elements relies on the
Mathematica packages FeynRules Alloul:2013bka , FeynArts Hahn:2000kx ,
FormCalc Hahn:1998yk ; Hahn:2016ebn , Package-X Patel:2015tea  and LiteRed
Lee:2013mka  as well as VSM.

The interference term between the tree-level matrix elements of $h(p_{h})\to
b(p_{1})\bar{b}(p_{2})g(p_{3})$ with an insertion of $Q_{bG}$ and the
corresponding SM contribution gives rise to a contribution of ${\cal
O}(y_{b}^{2}\hskip 0.7113pt\alpha_{s}^{2}\hskip 0.7113ptC_{bG})$. It takes the
following form

${\cal S}_{b\bar{b}g}^{(0)}=2\hskip 0.7113pt{\rm Re}\,\Big{\langle}{\cal
M}_{y_{b},b\bar{b}g}^{(0)}\,\Big{|}\,{\cal
M}_{C_{bG},b\bar{b}g}^{(0)}\Big{\rangle}=8\hskip 0.7113pt\alpha_{s}^{2}\hskip
0.7113ptC_{F}\hskip 0.7113ptC_{A}\hskip 0.7113ptm_{h}^{2}\hskip
1.42262pty_{b}^{2}\hskip 2.13394pt\frac{{\rm
Re}\left(C_{bG}\right)}{\Lambda^{2}}\,.$ (23)

Here ${\cal M}_{y_{b},b\bar{b}g}^{(0)}$ and ${\cal
M}_{C_{bG},b\bar{b}g}^{(0)}$ denote the SM and SMEFT tree-level $h\to
b\bar{b}g$ matrix element, respectively, and $C_{F}=4/3$ and $C_{A}=3$ are the
relevant colour factors.

At ${\cal O}(y_{b}^{2}\hskip 0.7113pt\alpha_{s}^{3}\hskip 0.7113ptC_{bG})$ the
interference term between the $Q_{bG}$ and the SM contribution to $h(p_{h})\to
b(p_{1})\bar{b}(p_{2})g(p_{3})$ can be written as follows

${\cal S}_{b\bar{b}g}^{(1)}=2\,{\rm Re}\,\Big{\langle}{\cal
M}_{y_{b},b\bar{b}g}^{(1)}\,\Big{|}\,{\cal
M}_{C_{bG},b\bar{b}g}^{(0)}\Big{\rangle}+2\,{\rm Re}\,\Big{\langle}{\cal
M}_{y_{b},b\bar{b}g}^{(0)}\,\Big{|}\,{\cal
M}_{C_{bG},b\bar{b}g}^{(1)}\Big{\rangle}\,,$ (24)

where ${\cal M}_{y_{b},b\bar{b}g}^{(1)}$ and ${\cal
M}_{C_{bG},b\bar{b}g}^{(1)}$ denote the SM and SMEFT one-loop $h\to b\bar{b}g$
matrix element. To write our result for ${\cal S}_{b\bar{b}g}^{(1)}$ in a
compact form we will use

$y_{ij}=\frac{2\hskip 0.7113ptp_{i}\cdot p_{j}}{m_{h}^{2}}\,,$ (25)

i.e. twice the dot-product of two external four-momenta $p_{i}$ and $p_{j}$
divided by the Higgs boson mass squared. For the sum of terms in (24) we
obtain

$\displaystyle\begin{split}{\cal
S}_{b\bar{b}g}^{(1)}=\frac{\alpha_{s}}{2\pi}\,N_{\epsilon}\,{\cal
S}_{b\bar{b}g}^{(0)}\,\Bigg{\\{}&\displaystyle-\frac{2\hskip
0.7113ptC_{F}+C_{A}}{\epsilon^{2}}-\frac{1}{\epsilon}\,\big{[}\left(C_{A}-2\hskip
0.7113ptC_{F}\right)\,L_{12}-C_{A}\left(L_{13}+L_{23}\right)+\gamma_{g}+2\hskip
0.35565pt\gamma_{q}\big{]}\\\\[5.69054pt] &\displaystyle+\left(C_{A}-2\hskip
0.7113ptC_{F}\right)\left[\frac{1}{2}\hskip 1.42262ptR\hskip
0.7113pt(y_{12},y_{13})+\frac{1}{2}\hskip 1.42262ptR\hskip
0.7113pt(y_{12},y_{23})+\frac{L_{12}^{2}}{2}-L_{12}\right]\\\\[5.69054pt]
&\displaystyle-C_{A}\left[\frac{1}{2}\hskip 1.42262ptR\hskip
0.7113pt(y_{13},y_{23})+\frac{L_{13}^{2}}{2}+\frac{L_{23}^{2}}{2}\right]-C_{F}\left(1-L_{13}-L_{23}\right)\\\\[5.69054pt]
&\displaystyle+3\hskip 1.42262pt\big{(}2\hskip
0.7113ptC_{F}+C_{A}\big{)}\hskip
1.42262pt\zeta_{2}+\frac{1}{2}\left(C_{A}-C_{F}\right)\left(y_{13}+y_{23}\right)+\gamma_{bG}+\left(\gamma_{g}+2\hskip
0.35565pt\gamma_{q}+\gamma_{bG}\right)L\hskip
1.42262pt\Bigg{\\}}\,.\end{split}$ (26)

Here we have defined

$N_{\epsilon}=\frac{e^{\epsilon\hskip 0.7113pt\gamma_{E}}}{\Gamma\hskip
0.7113pt(1-\epsilon)}\hskip
0.7113pt\left(\frac{\mu^{2}}{m_{h}^{2}}\right)^{\epsilon}=1+L\hskip
0.7113pt\epsilon+\left[L^{2}-\zeta_{2}\right]\frac{\epsilon^{2}}{2}\,,$ (27)

where $\Gamma\hskip 0.7113pt(z)$ is the gamma function, $\gamma_{E}\simeq
0.577$ is the Euler-Mascheroni constant, $\zeta_{2}=\pi^{2}/6$ is the Riemann
Zeta function of two and

$L_{ij}=\ln y_{ij}\,,\qquad L=\ln\frac{\mu^{2}}{m_{h}^{2}}\,.$ (28)

We have furthermore introduced

$R(x,y)={\rm Li}_{2}\hskip 0.7113pt(1-x)+{\rm Li}_{2}\hskip 0.7113pt(1-y)+\ln
x\,\ln y-\zeta_{2}\,,$ (29)

with ${\rm Li}_{2}(z)$ the dilogarithm. The anomalous dimensions $\gamma_{g}$
and $\gamma_{q}$ appearing in (26) can for example be found in the classic
work Catani:1998bh  on the IR structure of QCD scattering amplitudes, while
calculations of the anomalous dimension of the operator $Q_{bG}$ have been
performed in Misiak:1994zw ; Gorbahn:2005sa . The needed anomalous dimensions
read

$\gamma_{g}=\frac{11\hskip
0.7113ptC_{A}}{6}-\frac{N_{f}}{3}\,,\qquad\gamma_{q}=\frac{3\hskip
0.7113ptC_{F}}{2}\,,\qquad\gamma_{bG}=2\hskip 0.7113ptC_{A}-8\hskip
0.7113ptC_{F}\,,$ (30)

where $N_{f}$ denotes the number of active quark flavours. We add that as a
cross-check of our SMEFT calculation, we have also computed the interference
between the tree-level and one-loop $h\to b\bar{b}g$ matrix elements in the
SM. Our SM calculation reproduces the expression given in (A.6) of the
publication DelDuca:2015zqa .555The terms $\ln y_{ij}^{2}$ in (A.6) of
DelDuca:2015zqa  are obvious misprints that should in fact read
$\ln^{2}y_{ij}$.

We also need the ${\cal O}(y_{b}^{2}\hskip 0.7113pt\alpha_{s}^{3}\hskip
0.7113ptC_{bG})$ interferences for all relevant Higgs to four parton
scattering processes. To write these squared tree-level amplitudes in a
compact form we introduce

$y_{ijk}=y_{ij}+y_{ik}+y_{jk}\,.$ (31)

For $h(p_{h})\to b(p_{1})\bar{b}(p_{2})q(p_{3})\bar{q}(p_{4})$, we obtain

$\begin{split}{\cal S}_{b\bar{b}q\bar{q}}^{(0)}=2\hskip 0.7113pt{\rm
Re}\,\Big{\langle}{\cal M}_{y_{b},b\bar{b}q\bar{q}}^{(0)}\,\Big{|}\,{\cal
M}_{C_{bG},b\bar{b}q\bar{q}}^{(0)}\Big{\rangle}=2\hskip 0.7113pt\pi\hskip
0.7113pt\alpha_{s}\hskip 1.42262pt{\cal S}_{b\bar{b}g}^{(0)}\hskip
1.42262pt\frac{1}{m_{h}^{2}}\hskip
1.42262ptg_{b\bar{b}q\bar{q}}\left(p_{1},p_{2},p_{3},p_{4}\right)\,,\end{split}$
(32)

where

$\displaystyle\begin{split}g_{b\bar{b}q\bar{q}}\left(p_{1},p_{2},p_{3},p_{4}\right)&=-\frac{y_{13}}{y_{234}}+\frac{y_{13}}{y_{34}y_{134}}-\frac{2y_{24}y_{13}}{y_{134}y_{234}}-\frac{\left(y_{14}y_{23}+3y_{13}y_{24}\right)y_{13}}{y_{34}y_{134}y_{234}}\\\\[5.69054pt]
&-\frac{y_{24}\left(y_{13}y_{24}-y_{14}y_{23}\right)y_{13}}{y_{34}^{2}y_{134}y_{234}}+(1\leftrightarrow
2)+(3\leftrightarrow 4)+(1\leftrightarrow 2,3\leftrightarrow 4)\,.\hskip
22.76219pt\end{split}$ (33)

In the case of $h(p_{h})\to b(p_{1})\bar{b}(p_{2})b(p_{3})\bar{b}(p_{4})$, we
find

$\begin{split}{\cal S}_{b\bar{b}b\bar{b}}^{(0)}&=2\hskip 0.7113pt{\rm
Re}\,\Big{\langle}{\cal M}_{y_{b},b\bar{b}b\bar{b}}^{(0)}\,\Big{|}\,{\cal
M}_{C_{bG},b\bar{b}b\bar{b}}^{(0)}\Big{\rangle}\\\\[5.69054pt] &=2\hskip
0.7113pt\pi\hskip 0.7113pt\alpha_{s}\hskip 1.42262pt{\cal
S}_{b\bar{b}g}^{(0)}\hskip
1.42262pt\frac{1}{m_{h}^{2}}\,\bigg{[}\big{(}C_{A}-2\hskip
0.7113ptC_{F}\big{)}\hskip
1.42262ptf_{b\bar{b}b\bar{b}}\left(p_{1},p_{2},p_{3},p_{4}\right)+g_{b\bar{b}b\bar{b}}\left(p_{1},p_{2},p_{3},p_{4}\right)\bigg{]}\,,\end{split}$
(34)

where

$\displaystyle f_{b\bar{b}b\bar{b}}\left(p_{1},p_{2},p_{3},p_{4}\right)$
$\displaystyle=\frac{y_{12}^{5}}{y_{14}y_{123}y_{124}y_{134}y_{234}}+\frac{y_{12}^{5}}{y_{23}y_{123}y_{124}y_{134}y_{234}}+\frac{\left(4y_{13}+y_{23}+2y_{24}-1\right)y_{12}^{4}}{y_{14}y_{123}y_{124}y_{134}y_{234}}$
$\displaystyle+\frac{\left(y_{13}+y_{14}+4y_{24}-1\right)y_{12}^{4}}{y_{23}y_{123}y_{124}y_{134}y_{234}}+\frac{y_{13}\left(y_{23}+3y_{24}-4\right)y_{12}^{3}}{y_{14}y_{123}y_{124}y_{134}y_{234}}-\frac{\left(2y_{13}+y_{14}+y_{23}+y_{24}+2y_{34}\right)y_{12}^{3}}{y_{123}y_{124}y_{134}y_{234}}$
$\displaystyle+\frac{\left(y_{24}\left(2y_{13}+y_{14}+6y_{24}-4\right)-y_{14}y_{34}\right)y_{12}^{3}}{y_{23}y_{123}y_{124}y_{134}y_{234}}+\frac{\left(y_{24}\left(y_{23}+y_{24}-1\right)-y_{23}y_{34}\right)y_{12}^{3}}{y_{14}y_{123}y_{124}y_{134}y_{234}}$
$\displaystyle+\frac{6y_{13}^{2}y_{12}^{3}}{y_{14}y_{123}y_{124}y_{134}y_{234}}+\frac{y_{13}y_{24}^{2}y_{12}^{3}}{y_{23}y_{34}y_{123}y_{124}y_{134}y_{234}}+\frac{y_{13}^{2}y_{24}y_{12}^{3}}{y_{14}y_{34}y_{123}y_{124}y_{134}y_{234}}$
$\displaystyle-\frac{\left(y_{13}y_{23}-2y_{14}y_{23}+2y_{13}y_{24}+y_{14}y_{24}\right)y_{12}^{3}}{y_{34}y_{123}y_{124}y_{134}y_{234}}-\frac{y_{13}^{2}y_{12}^{2}}{y_{123}y_{124}y_{134}y_{234}}-\frac{y_{13}^{2}y_{24}y_{12}^{2}}{y_{23}y_{123}y_{124}y_{134}y_{234}}$
$\displaystyle+\frac{y_{13}^{2}\left(3y_{24}-4\right)y_{12}^{2}}{y_{14}y_{123}y_{124}y_{134}y_{234}}-\frac{y_{24}^{2}y_{34}y_{12}^{2}}{y_{14}y_{123}y_{124}y_{134}y_{234}}+\frac{\left(4y_{24}^{3}-4y_{24}^{2}-y_{34}\left(5y_{14}+y_{34}\right)y_{24}+2y_{14}y_{34}^{2}\right)y_{12}^{2}}{y_{23}y_{123}y_{124}y_{134}y_{234}}$
$\displaystyle+\frac{y_{13}\left(3y_{24}^{2}-4y_{34}y_{24}-y_{34}+y_{14}\left(y_{34}-5y_{24}\right)\right)y_{12}^{2}}{y_{23}y_{123}y_{124}y_{134}y_{234}}-\frac{y_{13}\left(3y_{14}+9y_{23}+4\left(3y_{24}+y_{34}\right)\right)y_{12}^{2}}{y_{123}y_{124}y_{134}y_{234}}$
$\displaystyle-\frac{y_{13}\left(5y_{23}\left(y_{24}+y_{34}\right)+y_{24}\left(y_{24}+5y_{34}\right)\right)y_{12}^{2}}{y_{14}y_{123}y_{124}y_{134}y_{234}}+\frac{\left(y_{14}^{2}-9y_{24}y_{14}+y_{23}^{2}-y_{24}\left(3y_{23}+y_{24}+5y_{34}\right)\right)y_{12}^{2}}{y_{123}y_{124}y_{134}y_{234}}$
$\displaystyle+\frac{3y_{13}^{3}y_{12}^{2}}{y_{14}y_{123}y_{124}y_{134}y_{234}}-\frac{2y_{13}^{2}y_{23}y_{12}^{2}}{y_{34}y_{123}y_{124}y_{134}y_{234}}+\frac{y_{13}y_{24}^{2}\left(3\left(y_{13}+y_{14}\right)+2y_{24}\right)y_{12}^{2}}{y_{23}y_{34}y_{123}y_{124}y_{134}y_{234}}$
$\displaystyle-\frac{y_{13}\left(3y_{23}^{2}+5\left(y_{14}+y_{23}\right)y_{24}\right)y_{12}^{2}}{y_{34}y_{123}y_{124}y_{134}y_{234}}+\frac{y_{14}\left(-2y_{24}^{2}-3y_{14}y_{24}+2y_{23}\left(y_{14}+y_{23}\right)\right)y_{12}^{2}}{y_{34}y_{123}y_{124}y_{134}y_{234}}$
$\displaystyle+\frac{y_{13}^{2}y_{24}\left(2y_{13}+3\left(y_{23}+y_{24}\right)\right)y_{12}^{2}}{y_{14}y_{34}y_{123}y_{124}y_{134}y_{234}}-\frac{y_{13}^{3}y_{12}}{y_{123}y_{124}y_{134}y_{234}}+\frac{\left(y_{24}-1\right)y_{24}^{3}y_{12}}{y_{23}y_{123}y_{124}y_{134}y_{234}}$
$\displaystyle+\frac{y_{13}y_{24}\left(y_{13}+y_{24}\right)\left(2y_{13}+y_{24}\right)y_{12}}{y_{14}y_{123}y_{124}y_{134}y_{234}}-\frac{y_{13}\left(y_{24}^{2}+4\left(y_{14}+y_{23}\right)y_{24}+6y_{14}y_{23}\right)y_{12}}{y_{123}y_{124}y_{134}y_{234}}$
$\displaystyle+\frac{3y_{13}y_{24}^{2}\left(y_{24}-y_{34}\right)y_{12}}{y_{23}y_{123}y_{124}y_{134}y_{234}}+\frac{3y_{13}^{2}y_{24}\left(y_{24}-y_{34}\right)y_{12}}{y_{23}y_{123}y_{124}y_{134}y_{234}}+\frac{y_{13}^{3}\left(y_{24}-y_{34}\right)y_{12}}{y_{23}y_{123}y_{124}y_{134}y_{234}}$
$\displaystyle-\frac{y_{13}^{2}\left(3y_{14}+4y_{23}+y_{24}+4y_{34}\right)y_{12}}{y_{123}y_{124}y_{134}y_{234}}-\frac{y_{24}\left(6y_{14}y_{23}+3y_{24}y_{23}+4y_{14}y_{24}+4y_{24}y_{34}\right)y_{12}}{y_{123}y_{124}y_{134}y_{234}}$
$\displaystyle-\frac{3y_{14}y_{23}y_{24}^{2}y_{12}}{y_{34}y_{123}y_{124}y_{134}y_{234}}+\frac{2y_{13}^{3}y_{24}^{2}y_{12}}{y_{23}y_{34}y_{123}y_{124}y_{134}y_{234}}-\frac{y_{13}^{3}\left(y_{23}-3y_{24}\right)y_{12}}{y_{34}y_{123}y_{124}y_{134}y_{234}}$
$\displaystyle+\frac{y_{13}^{4}y_{24}y_{12}}{y_{14}y_{34}y_{123}y_{124}y_{134}y_{234}}+\frac{y_{13}y_{24}^{3}\left(2y_{14}+y_{24}\right)y_{12}}{y_{23}y_{34}y_{123}y_{124}y_{134}y_{234}}+\frac{y_{13}^{2}y_{24}^{2}\left(3y_{23}+2y_{24}\right)y_{12}}{y_{14}y_{34}y_{123}y_{124}y_{134}y_{234}}$
$\displaystyle+\frac{y_{13}^{2}y_{24}^{2}\left(3y_{14}+5y_{24}\right)y_{12}}{y_{23}y_{34}y_{123}y_{124}y_{134}y_{234}}+\frac{y_{13}^{3}y_{24}\left(2y_{23}+5y_{24}\right)y_{12}}{y_{14}y_{34}y_{123}y_{124}y_{134}y_{234}}+\frac{y_{13}^{2}\left(6y_{24}^{2}-2y_{14}y_{24}-3y_{14}y_{23}\right)y_{12}}{y_{34}y_{123}y_{124}y_{134}y_{234}}$
$\displaystyle+\frac{y_{13}y_{24}\left(y_{24}\left(3y_{24}-2y_{23}\right)-8y_{14}y_{23}\right)y_{12}}{y_{34}y_{123}y_{124}y_{134}y_{234}}+\frac{2y_{13}y_{24}\left(y_{13}^{2}+y_{24}^{2}\right)}{y_{123}y_{124}y_{134}y_{234}}+\frac{y_{13}^{2}y_{24}^{3}\left(y_{13}+y_{24}\right)}{y_{23}y_{34}y_{123}y_{124}y_{134}y_{234}}$
$\displaystyle+\frac{y_{13}^{3}y_{24}^{2}\left(y_{13}+y_{24}\right)}{y_{14}y_{34}y_{123}y_{124}y_{134}y_{234}}+\frac{y_{13}y_{24}\left(y_{13}+y_{24}\right)\left(y_{13}^{2}+3y_{24}y_{13}+y_{24}^{2}\right)}{y_{34}y_{123}y_{124}y_{134}y_{234}}$
$\displaystyle+(1\leftrightarrow 3)+(2\leftrightarrow 4)+(1\leftrightarrow
3,2\leftrightarrow 4)\,.$ (35)

Furthermore,

$g_{b\bar{b}b\bar{b}}\left(p_{1},p_{2},p_{3},p_{4}\right)=g_{b\bar{b}q\bar{q}}\left(p_{1},p_{2},p_{3},p_{4}\right)+(1\leftrightarrow
3)+(2\leftrightarrow 4)+(1\leftrightarrow 3,2\leftrightarrow 4)\,,$ (36)

with $g_{b\bar{b}q\bar{q}}\left(p_{1},p_{2},p_{3},p_{4}\right)$ already
defined in (33).

For the process $h(p_{h})\to b(p_{1})\bar{b}(p_{2})g(p_{3})g(p_{4})$, we
finally get

$\begin{split}{\cal S}_{b\bar{b}gg}^{(0)}&=2\hskip 0.7113pt{\rm
Re}\,\Big{\langle}{\cal M}_{y_{b},b\bar{b}gg}^{(0)}\,\Big{|}\,{\cal
M}_{C_{bG},b\bar{b}gg}^{(0)}\Big{\rangle}\\\\[5.69054pt] &=2\hskip
0.7113pt\pi\hskip 0.7113pt\alpha_{s}\hskip 1.42262pt{\cal
S}_{b\bar{b}g}^{(0)}\hskip 1.42262pt\frac{1}{m_{h}^{2}}\hskip
1.42262pt\Big{[}\,C_{A}\hskip
0.7113ptf_{b\bar{b}gg}\left(p_{1},p_{2},p_{3},p_{4}\right)+C_{F}\hskip
0.7113ptg_{b\bar{b}gg}\left(p_{1},p_{2},p_{3},p_{4}\right)\Big{]}\,.\end{split}$
(37)

Here

$\displaystyle f_{b\bar{b}gg}\left(p_{1},p_{2},p_{3},p_{4}\right)$
$\displaystyle=\frac{4y_{24}^{2}}{y_{23}y_{34}y_{234}}+\frac{y_{13}^{2}y_{24}^{2}}{2y_{14}y_{23}y_{34}y_{134}y_{234}}+\frac{\left(10y_{13}+9y_{34}\right)y_{24}}{y_{14}y_{134}y_{234}}$
$\displaystyle+\frac{\left(4y_{24}+19y_{34}-4\right)y_{24}}{y_{23}y_{134}y_{234}}+\frac{\left(4y_{14}\left(y_{23}+2\right)+y_{13}\left(12y_{24}+7\right)\right)y_{24}}{2y_{34}y_{134}y_{234}}+\frac{2y_{13}\left(y_{13}y_{24}-y_{14}y_{23}\right)y_{24}}{y_{34}^{2}y_{134}y_{234}}$
$\displaystyle-\frac{2y_{34}^{2}}{y_{13}y_{14}y_{134}}+\frac{y_{34}^{2}\left(y_{23}+y_{24}+y_{34}-1\right)}{y_{13}y_{14}y_{23}y_{134}}+\frac{y_{34}^{2}\left(y_{23}+3y_{24}+2y_{34}\right)}{y_{13}y_{14}y_{134}y_{234}}$
$\displaystyle+\frac{18y_{34}+2y_{24}\left(4y_{13}+2y_{23}+2y_{24}+2y_{34}+19\right)-3}{2y_{134}y_{234}}+\frac{12y_{14}y_{24}+5y_{34}\left(2y_{24}+3y_{34}-1\right)}{y_{13}y_{134}y_{234}}$
$\displaystyle+\frac{2y_{14}y_{24}\left(2y_{24}+3y_{34}-1\right)+y_{34}\left(3y_{34}^{2}+\left(9y_{24}-2\right)y_{34}+4\left(y_{24}-1\right)y_{24}\right)}{y_{13}y_{23}y_{134}y_{234}}$
$\displaystyle+\frac{\left(y_{13}y_{24}\left(6y_{24}+8y_{34}-1\right)+y_{34}\left(6y_{24}^{2}+2\left(7y_{34}-2\right)y_{24}+y_{34}\left(5y_{34}-3\right)\right)\right)}{2y_{14}y_{23}y_{134}y_{234}}$
$\displaystyle+(1\leftrightarrow 2)+(3\leftrightarrow 4)+(1\leftrightarrow
2,3\leftrightarrow 4)\,,$ (38)

while

$\displaystyle g_{b\bar{b}gg}\left(p_{1},p_{2},p_{3},p_{4}\right)$
$\displaystyle=-\frac{2\left(y_{23}+y_{24}+y_{34}-1\right)y_{34}^{2}}{y_{13}y_{14}y_{23}y_{134}}+\frac{4y_{34}^{2}}{y_{13}y_{14}y_{134}}-\frac{2\left(y_{23}+3y_{24}+2y_{34}\right)y_{34}^{2}}{y_{13}y_{14}y_{134}y_{234}}$
$\displaystyle+\frac{2y_{24}\left(-4y_{24}+\left(y_{34}-14\right)y_{34}+4\right)}{y_{23}y_{134}y_{234}}+\frac{2y_{24}\left(y_{34}\left(y_{24}+2y_{34}-7\right)+y_{13}\left(y_{24}+2y_{34}-5\right)\right)}{y_{14}y_{134}y_{234}}$
$\displaystyle+\frac{2\left(3y_{24}^{2}+\left(4y_{13}+4y_{14}+3y_{23}+13y_{34}-14\right)y_{24}+y_{34}\left(3y_{34}-8\right)+2\right)}{y_{134}y_{234}}$
$\displaystyle-\frac{2\left(2y_{14}y_{24}\left(2y_{24}+3y_{34}-1\right)+y_{34}\left(3y_{34}^{2}+\left(9y_{24}-2\right)y_{34}+4\left(y_{24}-1\right)y_{24}\right)\right)}{y_{13}y_{23}y_{134}y_{234}}$
$\displaystyle+\frac{2\left(y_{14}y_{24}\left(y_{23}+y_{34}-8\right)+y_{34}\left(y_{24}^{2}+\left(2y_{23}+2y_{34}-9\right)y_{24}+\left(y_{34}-13\right)y_{34}+5\right)\right)}{y_{13}y_{134}y_{234}}$
$\displaystyle-\frac{2\left(y_{13}y_{24}\left(y_{24}+2y_{34}\right)+y_{34}\left(2y_{24}^{2}+\left(5y_{34}-1\right)y_{24}+y_{34}\left(2y_{34}-1\right)\right)\right)}{y_{14}y_{23}y_{134}y_{234}}$
$\displaystyle+(1\leftrightarrow 2)+(3\leftrightarrow 4)+(1\leftrightarrow
2,3\leftrightarrow 4)\,.$ (39)

As a cross-check of our SMEFT computation we have also calculated the squared
matrix elements for the tree-level processes $h\to b\bar{b}q\bar{q}$, $h\to
b\bar{b}b\bar{b}$ and $h\to b\bar{b}gg$ within the SM. Our findings agree with
(A.8), (A.9), (A.11), (A.13), (A.14) and (A.16) of DelDuca:2015zqa .

In the case of the $h(p_{h})\to b(p_{1})\bar{b}(p_{2})$ transition, the first
non-zero contribution to the interference with an insertion of $Q_{bG}$ and
the corresponding SM contribution arises at ${\cal O}(y_{b}^{2}\hskip
0.7113pt\alpha_{s}^{3}\hskip 0.7113ptC_{bG})$. We find

${\cal S}_{b\bar{b}}^{(2)}=2\hskip 0.7113pt{\rm Re}\,\Big{\langle}{\cal
M}_{y_{b},b\bar{b}}^{(0)}\,\Big{|}\,{\cal
M}_{C_{bG},b\bar{b}}^{(2)}\Big{\rangle}=-\frac{\alpha_{s}}{256\hskip
0.7113pt\pi^{3}}\,N_{\epsilon}^{2}\hskip 1.42262pt{\cal
S}_{b\bar{b}g}^{(0)}\hskip 1.42262ptm_{h}^{2}\hskip 1.42262pt\gamma_{bG}\hskip
1.42262pt\big{(}7+2\hskip 0.35565ptL\big{)}\,,$ (40)

where ${\cal M}_{C_{bG},b\bar{b}}^{(2)}$ denotes the relevant two-loop $h\to
b\bar{b}$ matrix element in the SMEFT and the expressions for ${\cal
S}_{b\bar{b}g}^{(0)}$, $N_{\epsilon}$ and $\gamma_{bG}$ have been given in
(23), (27) and (30), respectively.

##  Appendix B Numerical impact of $C_{HG}$ and $C_{3G}$

In this appendix we assess the numerical impact of the dimension-six operators
$Q_{HG}$ and $Q_{3G}$ on the $h\to b\bar{b}$ decay and the $pp\to Zh$
production process. Since we show below that these operators contribute to
less than a percent at the level of total rates once existing experimental
limits on the relevant Wilson coefficients are taken into account, the effects
of $Q_{HG}$ and $Q_{3G}$ have been neglected in the phenomenological analysis
carried out in the main part of this work.

The leading correction to the inclusive $h\to b\bar{b}$ decay rate
proportional to the Wilson coefficient $C_{HG}$ has been calculated for
massive bottom quarks in Gauld:2016kuu . Example diagrams are shown in Figure
6. Employing the operator basis (2) and working to leading power in $m_{b}$
one finds at NNLO in QCD the simple expression

$\Gamma(h\to b\bar{b})_{\rm SMEFT}^{{\rm
NNLO},HG}=\left(\frac{\alpha_{s}}{\pi}\right)^{2}\left[\frac{19}{3}-2\hskip
0.35565pt\zeta_{2}+\frac{1}{3}\ln^{2}\left(\frac{\bar{m}_{b}^{2}}{m_{h}^{2}}\right)\right]c_{HG}\hskip
1.42262pt\Gamma(h\to b\bar{b})_{\rm SM}^{\rm LO}\,,$ (41)

where the definition of the Wilson coefficient $c_{HG}$ can be found in (16).
A recent global fit Ellis:2020unq  to the SMEFT including 34 dimension-six
operators reports the following marginalised 95% CL bound

$c_{HG}\in[-0.09,0.06]\,,$ (42)

on the relevant Wilson coefficient (cf. Table 6 of the latter work). Using
this limit together with the input parameters given at the beginning of
Section 4 it follows from (41) and (42) that at 95% CL the $Q_{HG}$
contribution to the inclusive $h\to b\bar{b}$ decay rate lies within

$\frac{\Gamma(h\to b\bar{b})_{\rm SMEFT}^{{\rm NNLO},HG}}{\Gamma(h\to
b\bar{b})_{\rm SM}^{\rm LO}}\in[-2.7,1.7]\cdot 10^{-3}\,.$ (43)

This result indicates that SMEFT effects arising from $Q_{HG}$ are
phenomenologically irrelevant in the case of the fully differential $h\to
b\bar{b}$ decay rate. This justifies that we have neglected such corrections
in the main part of this article.

Figure 6:  Examples of SMEFT contributions to the inclusive decay rate of
$h\to b\bar{b}$ with an insertion of $Q_{HG}$ (black square). The left (right)
diagram represents a one-loop (tree-level) contribution to the $h\to b\bar{b}$
($h\to b\bar{b}g$) transition.

Operator insertions of $Q_{3G}$ induce tree-level contributions to $h\to
b\bar{b}gg$ and one-loop corrections to $h\to b\bar{b}g$. After interfering
these channels with the corresponding SM amplitudes and integrating over the
four- and three-particle phase space, respectively, the combination of these
two types of contributions leads to a N3LO correction to the inclusive $h\to
b\bar{b}$ decay rate. We write the sum of these real and virtual corrections
as

$\Gamma(h\to b\bar{b})_{\rm SMEFT}^{{{\rm N}^{3}{\rm LO}},3G}=N_{3G}^{\rm
dec}\left(\frac{\alpha_{s}}{\pi}\right)^{3}\frac{m_{h}^{2}}{v^{2}}\hskip
1.42262ptc_{3G}\hskip 1.42262pt\Gamma(h\to b\bar{b})_{\rm SM}^{\rm LO}\,,$
(44)

where we have defined

$c_{3G}=\frac{v^{2}}{\Lambda^{2}}\hskip 1.42262ptC_{3G}\,.$ (45)

The marginalised 95% CL limit on the relevant Wilson coefficient reads
Ellis:2020unq 666The large negative values of $c_{3G}$ found in the work
Ellis:2020unq  can be traced back to the discrepancy between the measured
$t\bar{t}$ differential cross section CMS:2018htd  and the state-of-the-art SM
prediction Czakon:2017wor  at large values of the top-antitop invariant mass
$m_{t\bar{t}}$. Analyses of multijet data Krauss:2016ely ; Hirschi:2018etq ;
Goldouzian:2020wdq  suggest bounds of $|c_{3G}|\lesssim 0.2$.

$c_{3G}\in[-12.5,-4.1]\,.$ (46)

Plugging (46) into (44) then leads to

$\frac{\Gamma(h\to b\bar{b})_{\rm SMEFT}^{{{\rm N}^{3}{\rm
LO}},3G}}{\Gamma(h\to b\bar{b})_{\rm SM}^{\rm LO}}\in[-0.15,-0.05]\cdot
10^{-3}\hskip 0.7113ptN_{3G}^{\rm dec}\,.$ (47)

By performing an explicit calculation of the tree-level contributions to $h\to
b\bar{b}gg$ and the one-loop corrections to $h\to b\bar{b}g$ associated to
insertions of $Q_{3G}$, we find for the unknown constant $N_{3G}^{\rm dec}$
introduced in (44) the numerial value $N_{3G}^{\rm dec}=2.23$. This implies
that the relative corrections associated to the operator $Q_{3G}$ do not even
reach the level of a permille. Neglecting these corrections as done in our
SMEFT analysis of the fully differential $h\to b\bar{b}$ decay rate is
therefore fully justified from a phenomenological point of view.

The leading corrections to $Zh$ production associated to the operator $Q_{HG}$
result from the Feynman graphs displayed in Figure 7. These types of diagrams
have been calculated in the context of the SM in Brein:2011vx  working in the
limit of infinite top-quark mass. Using the results of the latter work one can
write the ${\cal O}(\alpha_{s}^{2}\hskip 1.42262ptC_{HG})$ corrections to the
inclusive $Zh$ production cross section in the following way

$\sigma(pp\to Zh)_{\rm SMEFT}^{{\rm
NNLO},HG}=3\left(\frac{\alpha_{s}}{\pi}\right)^{2}\hskip 1.42262pt\delta\hskip
1.42262ptc_{HG}\,\sigma(pp\to Zh)_{\rm SM}^{\rm LO}\,,$ (48)

where $\sigma(pp\to Zh)_{\rm SM}^{\rm LO}$ is the LO cross section and
$\delta$ encodes the sum of the QCD corrections denoted by $V_{\rm I}$ and
$R_{\rm I}$ in Brein:2011vx  with a factor of $(\alpha_{s}/\pi)^{2}$ stripped
off. We add that we have calculated the relevant radiative corrections that
give rise to $V_{\rm I}$ and $R_{\rm I}$ finding agreement with the latter
publication (see also Brein:2012ne ; Harlander:2018yio ). Notice that an
expression analogous to (48) also holds for the differential cross section.
From Figure 6 of the paper Brein:2011vx  one finds that at the LHC one has
$\delta=10.7$ for the measured mass of the Higgs boson. This corresponds to a
$1.4\%$ correction in the SM. Using the limit (42) in (48) we obtain

$\frac{\sigma(pp\to Zh)_{\rm SMEFT}^{{{\rm NNLO}},HG}}{\sigma(pp\to Zh)_{\rm
SM}^{\rm LO}}\in[-3.9,2.4]\cdot 10^{-3}\,,$ (49)

at 95% CL. This numerical result shows that it is an excellent approximation
to neglect contributions due to $Q_{HG}$ in the calculation of $Zh$ production
observables.

Figure 7:  Examples of SMEFT contributions to $Zh$ production involving an
insertion of $Q_{HG}$ (black square). The left and middle (right) diagram
represent(s) a one-loop (tree-level) contribution to the $q\bar{q}\to Zh$
($q\bar{q}\to Zhg$) process.

Insertions of the operator $Q_{3G}$ induce tree-level contributions to
$q\bar{q}\to Zhgg$ and one-loop corrections to $q\bar{q}\to Zhg$. Interfering
the SMEFT with the corresponding SM amplitudes and integrating over the
relevant phase spaces, one obtains a N3LO correction to the inclusive $Zh$
production cross section. Using naive dimensional analysis, we write this
correction as follows

$\sigma(pp\to Zh)_{\rm SMEFT}^{{{\rm N}^{3}{\rm LO}},3G}=N_{3G}^{\rm
prod}\left(\frac{\alpha_{s}}{\pi}\right)^{3}\hskip
1.42262ptc_{3G}\,\sigma(pp\to Zh)_{\rm SM}^{\rm LO}\,,$ (50)

where $N_{3G}^{\rm prod}$ is an undetermined numerical factor which is naively
of ${\cal O}(1)$. From (46) and (50) it follows that

$\frac{\sigma(pp\to Zh)_{\rm SMEFT}^{{{\rm N}^{3}{\rm LO}},3G}}{\sigma(pp\to
Zh)_{\rm SM}^{\rm LO}}\in[-0.58,-0.19]\cdot 10^{-3}\hskip 0.7113ptN_{3G}^{\rm
prod}\,.$ (51)

For $N_{3G}^{\rm prod}={\cal O}(10)$ this 95% CL bound implies a relative
correction of a few permille to the inclusive $pp\to Zh$ cross section. In
practice, the ${{\rm N}^{3}{\rm LO}}$ corrections (50) can therefore be
neglected.

## References

  * (1) ATLAS collaboration, _Observation of $H\rightarrow b\bar{b}$ decays and $VH$ production with the ATLAS detector_, _Phys. Lett. B_ 786 (2018) 59 [1808.08238].
  * (2) CMS collaboration, _Observation of Higgs boson decay to bottom quarks_ , _Phys. Rev. Lett._ 121 (2018) 121801 [1808.08242].
  * (3) ATLAS collaboration, _Projections for measurements of Higgs boson cross sections, branching ratios, coupling parameters and mass with the ATLAS detector at the HL-LHC_ , ATL-PHYS-PUB-2018-054.
  * (4) CMS collaboration, _Sensitivity projections for Higgs boson properties measurements at the HL-LHC_ , CMS-PAS-FTR-18-011.
  * (5) W. Buchmüller and D. Wyler, _Effective Lagrangian Analysis of New Interactions and Flavor Conservation_ , _Nucl. Phys._ B268 (1986) 621.
  * (6) B. Grzadkowski, M. Iskrzynski, M. Misiak and J. Rosiek, _Dimension-Six Terms in the Standard Model Lagrangian_ , _JHEP_ 10 (2010) 085 [1008.4884].
  * (7) I. Brivio and M. Trott, _The Standard Model as an Effective Field Theory_ , _Phys. Rept._ 793 (2019) 1 [1706.08945].
  * (8) K. Mimasu, V. Sanz and C. Williams, _Higher Order QCD predictions for Associated Higgs production with anomalous couplings to gauge bosons_ , _JHEP_ 08 (2016) 039 [1512.02572].
  * (9) C. Degrande, B. Fuks, K. Mawatari, K. Mimasu and V. Sanz, _Electroweak Higgs boson production in the standard model effective field theory beyond leading order in QCD_ , _Eur. Phys. J. C_ 77 (2017) 262 [1609.04833].
  * (10) R. Gauld, B.D. Pecjak and D.J. Scott, _One-loop corrections to $h\to b\bar{b}$ and $h\to\tau\bar{\tau}$ decays in the Standard Model Dimension-6 EFT: four-fermion operators and the large-$m_{t}$ limit_, _JHEP_ 05 (2016) 080 [1512.02508].
  * (11) R. Gauld, B.D. Pecjak and D.J. Scott, _QCD radiative corrections for $h\to b\bar{b}$ in the Standard Model Dimension-6 EFT_, _Phys. Rev. D_ 94 (2016) 074045 [1607.06354].
  * (12) J.M. Cullen, B.D. Pecjak and D.J. Scott, _NLO corrections to $h\to b\bar{b}$ decay in SMEFT_, _JHEP_ 08 (2019) 173 [1904.06358].
  * (13) J.M. Cullen and B.D. Pecjak, _Higgs decay to fermion pairs at NLO in SMEFT_ , _JHEP_ 11 (2020) 079 [2007.15238].
  * (14) W. Bizoń, F. Caola, K. Melnikov and R. Röntsch, _Anomalous couplings in associated VH production with Higgs boson decay to massive b quarks at NNLO in QCD_ , _Phys. Rev. D_ 105 (2022) 014023 [2106.06328].
  * (15) W. Astill, W. Bizoń, E. Re and G. Zanderighi, _NNLOPS accurate associated HZ production with $H\to b\overline{b}$ decay at NLO_, _JHEP_ 11 (2018) 157 [1804.08141].
  * (16) S. Alioli, A. Broggio, S. Kallweit, M.A. Lim and L. Rottoli, _Higgsstrahlung at NNLL ${}^{\prime}+$NNLO matched to parton showers in GENEVA_, _Phys. Rev. D_ 100 (2019) 096016 [1909.02026].
  * (17) W. Bizoń, E. Re and G. Zanderighi, _NNLOPS description of the $H\to b\overline{b}$ decay with MiNLO_, _JHEP_ 06 (2020) 006 [1912.09982].
  * (18) S. Zanoli, M. Chiesa, E. Re, M. Wiesemann and G. Zanderighi, _Next-to-next-to-leading order event generation for $VH$ production with $H\to b\bar{b}$ decay_, [2112.04168].
  * (19) F. Granata, J.M. Lindert, C. Oleari and S. Pozzorini, _NLO QCD+EW predictions for HV and HV+jet production including parton-shower effects_ , _JHEP_ 09 (2017) 012 [1706.03522].
  * (20) G.F. Giudice, C. Grojean, A. Pomarol and R. Rattazzi, _The Strongly-Interacting Light Higgs_ , _JHEP_ 06 (2007) 045 [hep-ph/0703164].
  * (21) J. Elias-Miro, J.R. Espinosa, E. Masso and A. Pomarol, _Higgs windows to new physics through d=6 operators: constraints and one-loop anomalous dimensions_ , _JHEP_ 11 (2013) 066 [1308.1879].
  * (22) G. D’Ambrosio, G.F. Giudice, G. Isidori and A. Strumia, _Minimal flavor violation: An Effective field theory approach_ , _Nucl. Phys._ B645 (2002) 155 [hep-ph/0207036].
  * (23) C. Anastasiou, F. Herzog and A. Lazopoulos, _The fully differential decay rate of a Higgs boson to bottom-quarks at NNLO in QCD_ , _JHEP_ 03 (2012) 035 [1110.2368].
  * (24) V. Del Duca, C. Duhr, G. Somogyi, F. Tramontano and Z. Trócsányi, _Higgs boson decay into b-quarks at NNLO accuracy_ , _JHEP_ 04 (2015) 036 [1501.07226].
  * (25) F. Caola, K. Melnikov and R. Röntsch, _Analytic results for decays of color singlets to $gg$ and $q\bar{q}$ final states at NNLO QCD with the nested soft-collinear subtraction scheme_, _Eur. Phys. J. C_ 79 (2019) 1013 [1907.05398].
  * (26) S.G. Gorishnii, A.L. Kataev and S.A. Larin, _The Width of Higgs Boson Decay Into Hadrons: Three Loop Corrections of Strong Interactions_ , _Sov. J. Nucl. Phys._ 40 (1984) 329.
  * (27) S.G. Gorishnii, A.L. Kataev, S.A. Larin and L.R. Surguladze, _Corrected Three Loop QCD Correction to the Correlator of the Quark Scalar Currents and $\Gamma_{\rm tot}(H^{0}\to{\rm hadrons})$_, _Mod. Phys. Lett. A_ 5 (1990) 2703.
  * (28) S.G. Gorishnii, A.L. Kataev, S.A. Larin and L.R. Surguladze, _Scheme dependence of the next to next-to-leading QCD corrections to $\Gamma_{\rm tot}(H^{0}\to{\rm hadrons})$ and the spurious QCD infrared fixed point_, _Phys. Rev. D_ 43 (1991) 1633.
  * (29) K.G. Chetyrkin, _Correlator of the quark scalar currents and $\Gamma_{\rm tot}(H\to{\rm hadrons})$ at ${\cal O}(\alpha_{s}^{3})$ in pQCD_, _Phys. Lett. B_ 390 (1997) 309 [hep-ph/9608318].
  * (30) P.A. Baikov, K.G. Chetyrkin and J.H. Kühn, _Scalar correlator at ${\cal O}(\alpha_{s}^{4})$, Higgs decay into $b$-quarks and bounds on the light quark masses_, _Phys. Rev. Lett._ 96 (2006) 012003 [hep-ph/0511063].
  * (31) F. Herzog, B. Ruijl, T. Ueda, J.A.M. Vermaseren and A. Vogt, _On Higgs decays to hadrons and the R-ratio at N 4LO_, _JHEP_ 08 (2017) 113 [1707.01044].
  * (32) R. Contino, M. Ghezzi, C. Grojean, M. Mühlleitner and M. Spira, _Effective Lagrangian for a light Higgs-like scalar_ , _JHEP_ 07 (2013) 035 [1303.3876].
  * (33) R. Contino, M. Ghezzi, C. Grojean, M. Mühlleitner and M. Spira, _eHDECAY: an Implementation of the Higgs Effective Lagrangian into HDECAY_ , _Comput. Phys. Commun._ 185 (2014) 3412 [1403.3381].
  * (34) U. Haisch and G. Koole, _Beautiful and charming chromodipole moments_ , _JHEP_ 09 (2021) 133 [2106.01289].
  * (35) A. Hayreter and G. Valencia, _Constraints on anomalous color dipole operators from Higgs boson production at the LHC_ , _Phys. Rev. D_ 88 (2013) 034033 [1304.6976].
  * (36) J. Bramante, A. Delgado, L. Lehman and A. Martin, _Boosted Higgses from chromomagnetic $b$’s: $b\bar{b}h$ at high luminosity_, _Phys. Rev. D_ 93 (2016) 053001 [1410.3484].
  * (37) K.G. Chetyrkin and A. Kwiatkowski, _Second order QCD corrections to scalar and pseudoscalar Higgs decays into massive bottom quarks_ , _Nucl. Phys. B_ 461 (1996) 3 [hep-ph/9505358].
  * (38) S.A. Larin, T. van Ritbergen and J.A.M. Vermaseren, _The Large top quark mass expansion for Higgs boson decays into bottom quarks and into gluons_ , _Phys. Lett. B_ 362 (1995) 134 [hep-ph/9506465].
  * (39) W. Bernreuther, L. Chen and Z.-G. Si, _Differential decay rates of CP-even and CP-odd Higgs bosons to top and bottom quarks at NNLO QCD_ , _JHEP_ 07 (2018) 159 [1805.06658].
  * (40) A. Behring and W. Bizoń, _Higgs decay into massive b-quarks at NNLO QCD in the nested soft-collinear subtraction scheme_ , _JHEP_ 01 (2020) 189 [1911.11524].
  * (41) A. Behring, W. Bizoń, F. Caola, K. Melnikov and R. Röntsch, _Bottom quark mass effects in associated $WH$ production with the $H\to b\bar{b}$ decay through NNLO QCD_, _Phys. Rev. D_ 101 (2020) 114012 [2003.08321].
  * (42) G. Somogyi and F. Tramontano, _Fully exclusive heavy quark-antiquark pair production from a colourless initial state at NNLO in QCD_ , _JHEP_ 11 (2020) 142 [2007.15015].
  * (43) R. Akhoury, H. Wang and O.I. Yakovlev, _On the Resummation of large QCD logarithms in $H\to\gamma\gamma$ decay_, _Phys. Rev. D_ 64 (2001) 113008 [hep-ph/0102105].
  * (44) T. Liu and A.A. Penin, _High-Energy Limit of QCD beyond the Sudakov Approximation_ , _Phys. Rev. Lett._ 119 (2017) 262001 [1709.01092].
  * (45) Z.L. Liu and M. Neubert, _Factorization at subleading power and endpoint-divergent convolutions in $h\to\gamma\gamma$ decay_, _JHEP_ 04 (2020) 033 [1912.08818].
  * (46) J. Wang, _Resummation of double logarithms in loop-induced processes with effective field theory_ , [1912.09920].
  * (47) O. Brein, R. Harlander, M. Wiesemann and T. Zirke, _Top-Quark Mediated Effects in Hadronic Higgs-Strahlung_ , _Eur. Phys. J. C_ 72 (2012) 1868 [1111.0761].
  * (48) K. Hamilton, P. Nason and G. Zanderighi, _MINLO: Multi-Scale Improved NLO_ , _JHEP_ 10 (2012) 155 [1206.3572].
  * (49) K. Hamilton, P. Nason, C. Oleari and G. Zanderighi, _Merging H/W/Z + 0 and 1 jet at NLO with no merging scale: a path to parton shower + NNLO matching_ , _JHEP_ 05 (2013) 082 [1212.4504].
  * (50) P. Nason, _A New method for combining NLO QCD with shower Monte Carlo algorithms_ , _JHEP_ 11 (2004) 040 [hep-ph/0409146].
  * (51) S. Frixione, P. Nason and C. Oleari, _Matching NLO QCD computations with Parton Shower simulations: the POWHEG method_ , _JHEP_ 11 (2007) 070 [0709.2092].
  * (52) S. Frixione, Z. Kunszt and A. Signer, _Three jet cross-sections to next-to-leading order_ , _Nucl. Phys. B_ 467 (1996) 399 [hep-ph/9512328].
  * (53) S. Frixione, _A General approach to jet cross-sections in QCD_ , _Nucl. Phys. B_ 507 (1997) 295 [hep-ph/9706545].
  * (54) S. Alioli, P. Nason, C. Oleari and E. Re, _A general framework for implementing NLO calculations in shower Monte Carlo programs: the POWHEG BOX_ , _JHEP_ 06 (2010) 043 [1002.2581].
  * (55) Particle Data Group collaboration, _Review of Particle Physics_ , _PTEP_ 2020 (2020) 083C01.
  * (56) R.D. Ball, V. Bertone, S. Carrazza et al., _Parton distributions from high-precision collider data_ , _Eur. Phys. J. C_ 77 (2017) 663 [1706.00428].
  * (57) T. Sjöstrand, S. Ask, J.R. Christiansen, R. Corke, N. Desai, P. Ilten et al., _An introduction to PYTHIA 8.2_ , _Comput. Phys. Commun._ 191 (2015) 159 [1410.3012].
  * (58) P. Skands, S. Carrazza and J. Rojo, _Tuning PYTHIA 8.1: the Monash 2013 Tune_ , _Eur. Phys. J. C_ 74 (2014) 3024 [1404.5630].
  * (59) J. Ellis, M. Madigan, K. Mimasu, V. Sanz and T. You, _Top, Higgs, Diboson and Electroweak Fit to the Standard Model Effective Field Theory_ , _JHEP_ 04 (2021) 279 [2012.02779].
  * (60) ATLAS collaboration, _Measurements of the production cross-section for a $Z$ boson in association with $b$-jets in proton-proton collisions at $\sqrt{s}=13\,{\rm TeV}$ with the ATLAS detector_, _JHEP_ 07 (2020) 044 [2003.11960].
  * (61) ATLAS collaboration, _Search for new resonances in mass distributions of jet pairs using $139\,{\rm fb}^{-1}$ of $pp$ collisions at $\sqrt{s}=13\,{\rm TeV}$ with the ATLAS detector_, _JHEP_ 03 (2020) 145 [1910.08447].
  * (62) M. Cacciari, G.P. Salam and G. Soyez, _The anti- $k_{t}$ jet clustering algorithm_, _JHEP_ 04 (2008) 063 [0802.1189].
  * (63) M. Cacciari, G.P. Salam and G. Soyez, _FastJet User Manual_ , _Eur. Phys. J. C_ 72 (2012) 1896 [1111.6097].
  * (64) J.M. Butterworth, A.R. Davison, M. Rubin and G.P. Salam, _Jet substructure as a new Higgs search channel at the LHC_ , _Phys. Rev. Lett._ 100 (2008) 242001 [0802.2470].
  * (65) J.R. Andersen et al., _Les Houches 2015: Physics at TeV Colliders Standard Model Working Group Report_ , in _9th Les Houches Workshop on Physics at TeV Colliders_ , 5, 2016 [1605.04692].
  * (66) N. Berger et al., _Simplified Template Cross Sections - Stage 1.1_ , [1906.02754].
  * (67) S. Amoroso et al., _Les Houches 2019: Physics at TeV Colliders: Standard Model Working Group Report_ , in _11th Les Houches Workshop on Physics at TeV Colliders: PhysTeV Les Houches_ , 3, 2020 [2003.01700].
  * (68) ATLAS collaboration, _Measurement of $VH$, $H\to b\bar{b}$ production as a function of the vector-boson transverse momentum in $13\,{\rm TeV}$ $pp$ collisions with the ATLAS detector_, _JHEP_ 05 (2019) 141 [1903.04618].
  * (69) ATLAS collaboration, _Measurements of $WH$ and $ZH$ production in the $H\rightarrow b\bar{b}$ decay channel in $pp$ collisions at $13\,{\rm TeV}$ with the ATLAS detector_, _Eur. Phys. J. C_ 81 (2021) 178 [2007.02873].
  * (70) P.F. Monni, P. Nason, E. Re, M. Wiesemann and G. Zanderighi, _MiNNLO PS: a new method to match NNLO QCD to parton showers_, _JHEP_ 05 (2020) 143 [1908.06987].
  * (71) P.F. Monni, E. Re and M. Wiesemann, _MiNNLO ${}_{\text{PS}}$: optimizing $2\rightarrow 1$ hadronic processes_, _Eur. Phys. J. C_ 80 (2020) 1075 [2006.04133].
  * (72) G. Coloretti, A. Gehrmann-De Ridder and C.T. Preuss, _QCD Predictions for Event-Shape Distributions in Hadronic Higgs Decays_ , [2202.07333].
  * (73) _The POWHEG BOX_ , http://powhegbox.mib.infn.it.
  * (74) A. Alloul, N.D. Christensen, C. Degrande, C. Duhr and B. Fuks, _FeynRules 2.0 - A complete toolbox for tree-level phenomenology_ , _Comput. Phys. Commun._ 185 (2014) 2250 [1310.1921].
  * (75) T. Hahn, _Generating Feynman diagrams and amplitudes with FeynArts 3_ , _Comput. Phys. Commun._ 140 (2001) 418 [hep-ph/0012260].
  * (76) T. Hahn and M. Perez-Victoria, _Automatized one loop calculations in four-dimensions and D-dimensions_ , _Comput. Phys. Commun._ 118 (1999) 153 [hep-ph/9807565].
  * (77) T. Hahn, S. Paßehr and C. Schappacher, _FormCalc 9 and Extensions_ , _PoS_ LL2016 (2016) 068 [1604.04611].
  * (78) H.H. Patel, _Package-X: A Mathematica package for the analytic calculation of one-loop integrals_ , _Comput. Phys. Commun._ 197 (2015) 276 [1503.01469].
  * (79) R.N. Lee, _LiteRed 1.4: a powerful tool for reduction of multiloop integrals_ , _J. Phys. Conf. Ser._ 523 (2014) 012059 [1310.1145].
  * (80) S. Catani, _The Singular behavior of QCD amplitudes at two loop order_ , _Phys. Lett. B_ 427 (1998) 161 [hep-ph/9802439].
  * (81) M. Misiak and M. Münz, _Two loop mixing of dimension five flavor changing operators_ , _Phys. Lett. B_ 344 (1995) 308 [hep-ph/9409454].
  * (82) M. Gorbahn, U. Haisch and M. Misiak, _Three-loop mixing of dipole operators_ , _Phys. Rev. Lett._ 95 (2005) 102004 [hep-ph/0504194].
  * (83) CMS collaboration, _Measurement of differential cross sections for the production of top quark pairs and of additional jets in lepton+jets events from $pp$ collisions at $\sqrt{s}=13\,{\rm TeV}$_, _Phys. Rev. D_ 97 (2018) 112003 [1803.08856].
  * (84) M. Czakon, D. Heymes, A. Mitov, D. Pagani, I. Tsinikos and M. Zaro, _Top-pair production at the LHC through NNLO QCD and NLO EW_ , _JHEP_ 10 (2017) 186 [1705.04105].
  * (85) F. Krauss, S. Kuttimalai and T. Plehn, _LHC multijet events as a probe for anomalous dimension-six gluon interactions_ , _Phys. Rev. D_ 95 (2017) 035024 [1611.00767].
  * (86) V. Hirschi, F. Maltoni, I. Tsinikos and E. Vryonidou, _Constraining anomalous gluon self-interactions at the LHC: a reappraisal_ , _JHEP_ 07 (2018) 093 [1806.04696].
  * (87) R. Goldouzian and M.D. Hildreth, _LHC dijet angular distributions as a probe for the dimension-six triple gluon vertex_ , _Phys. Lett. B_ 811 (2020) 135889 [2001.02736].
  * (88) O. Brein, R.V. Harlander and T.J.E. Zirke, _vh@nnlo - Higgs Strahlung at hadron colliders_ , _Comput. Phys. Commun._ 184 (2013) 998 [1210.5347].
  * (89) R.V. Harlander, J. Klappert, S. Liebler and L. Simon, _vh@nnlo-v2: New physics in Higgs Strahlung_ , _JHEP_ 05 (2018) 089 [1802.04817].